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BIDIRECTIONAL PLANT CANOPY REFLECTION MODELS 


DERIVED FROM THE RADIATION TRANSFER EQUATION 

By Donald R. Beeth 
Lyndon B. Johnson Space Center 

SUMMARY 


The basic approach to canopy modeling adopted in this paper is to solve the 
radiation transfer equation under assumptions appropriate for a canopy and thus 
to predict canopy- reflected radiation under arbitrary conditions of illumination . 
Consideration of the radiation transfer equation and its formal solution along a path 
shows that the basic problem of canopy modeling has two parts . The physical part 
of the problem is to define the phase function by determining the scattering behavior 
of a volume element that contains a collection of scattering elements . The mathemati- 
cal part of the problem is to solve the radiation transfer equation for the medium as 
a whole , not merely along a simple path as is done in the formal solution . Solving 
the problem for the medium as a whole requires the solution of an infinite set of inte- 
grodifferential equations . The phase function is defined by developing a technique 
to sum the scattering behavior of a collection of planar scattering elements . Coordi- 
nate transformations are developed to enable proper treatment of the potential bi- 
directional scattering behavior of the scattering elements. Canopy geometry is ex- 
pressed in terms of probability density functions over azimuth and zenith angles . 

Two approaches are presented for solving the mathematical problem . One 
approach is an extension of a Chandrasekhar technique . The approach concludes 
with eight integrodifferential equations that are reducible to six by transmission 
reciprocity. Because no programable solution is presented, an alternate approach 
is emphasized . The basic assumption for this second approach is that the field can 
be partitioned into a uniformly diffuse upwelling component , a uniformly diffuse 
downwelling component, and an attenuated specular component. This approach re- 
duces the radiation transfer equation to three easily solvable coupled, first-order 
linear differential equations of the Shuster or Kubelka-Munk type . The solution 
to these three equations gives an approximation to the field within the canopy and 
is used in the formal solution of the radiation transfer equation to obtain canopy bi- 
directional reflection. 


INTRODUCTION 


The objective for developing a reflective canopy model is to predict the up- 
welling monochromatic intensity from a canopy under specified conditions of illumi- 
nation , preferably in terms of physically measurable characteristics of the csinopy . 



These characteristics are the geometrical and optical properties of the components 
that make up a canopy . If these canopy component optical and geometric properties 
can be related to canopy type, stress, maturity, and yield, the model can be used 
to relate canopy-reflected radiation to these characteristics . This relationship would 
imply potential application to the acquisition and analysis of signature statistical 
properties . 

The approach presented in this paper results in three canopy models . The 
models vary in form and approach depending primarily on the combination of simpli- 
fying assumptions made . All models assume horizontally homogeneous canopies . 

In Model I and Model II , it is assumed that a first-order approximation of the field 
within the canopy is obtainable from modified Kubelka-Munk equations . In Model I , 
the acceptance of non-Lambertian leaves is possible if the non- Lambertian character 
of the leaves is known . Further , an azimuthal anisotropic distribution of leaves 
will be accepted if this anisotropic distribution is known . In Model II , the further 
simplifying assumption that the leaves are Lambertian scatterers is made . This 
assumption results in considerable simplification of the approach used and the calcu- 
lated expressions obtained. From the further assumption in Model II that the two 
sides of the leaves are optically indistinguishable, an interesting result is obtained: 
that all geometric factors affecting scattering of the diffuse field are multiplied by 
the difference between radiation reflected and radiation transmitted. This result 
would imply that changes in canopy geometry would have little effect in regions of 
canopy optical thickness sufficient to ensure that soil background has little effect. 

For example, wind-induced data scattering in these regions would be minimal. 

In Model III , which is the most general model , non-Lambertian leaves and 
azimuthal anisotropy are accepted and the Kubelka-Munk approximation is not re- 
quired . Although Model III is exact in principle , difficulties in programing make 
it , for the most part , a theoretical curiosity . However , the effect of such assump- 
tions as Lambertian leaves , azimuthal isotropy , and some well-defined zenith distri- 
bution density functions is not explored in this paper . Perhaps Model III could be 
made quite workable for some special situations , but this possibility is left for 
further investigation. 


SYMBOLS 


Ay Kubelka-Munk equation expansion coefficients 


A^^ expansion coefficients for expressing in terms of 

eigenfunction F^ as defined in equation (77) 


^ ^11 ^22 
a,b ,c ,d 


2 


integration frame limits 

Kubelka-Munk equation expansion coefficients 



B 

expression containing canopy geometric characteristics 
for canopies with Lambertian leaves and optically indis 
tinguishable sides 

b = ai2 = 


C,V,R 

transformation matrices 

^ ^ ^13 


~ ^23 


II 

Q 

derivative operator 

D’ 

1 2 

average value of — G (0^^) 

D(9n»'^) 

azimuthal density function 

dA 

element of horizontal area 

da 

magnitude of da 

da 

vector surface element 

<da> 

average projected area of a collection of surface elements 

da 

cross section of cylindrical volume element 

dan = da 


dI(C,k) 

differential change in intensity along 8 in the direction 
of k 

ds 

differential vector area on the surface of volume element 
AV 

dV 

volume element 

dn' d(f>’ = -dri' 


da' 

differential solid angle 

da. 

1 

differential solid angle containing incident radiation 

da 

s 

differential solid angle containing scattered radiation 
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E^(0) 


E^Jd.z) 


E^^<-d,z) 


E^^(s.z) 


irradiance of the specular beam measured perpendicular to 
the beam 

horizontal irradiance resulting from reduced specular flux 

irradiance on a horizontal plane at the upper surface of a 
canopy 

upwelling irradiance at level z in the a layer 


downwelling irradiance at level z in the a layer 


irradiance at level z resulting from the reduced incident 
specular flux 


E 


+ 


E 


E 


+ 


= %L^ 


F 


F 

a 




F 


($(t) 

G = ^(1 - cos 0) 
1 

o o 2 

g = (a^ - b^) 

±h 

I 

Kk) 

l(fi,k) 


upwelling irradiance 
downwelling irradiance 


fraction of beam obscured while passing through AV 
Shuster or Kubelka-Munk equation eigenfunction 

magnitude of incident intensity 

variable defined in equation (274) 
an arbitrary function defined in equation (174) 
variable defined in equations (99) and (100) 


integration over upper or lower hemispheres 

intensity where the functional dependence is suppressed 
upwelling intensity at the upper surface of a canopy 

A 

intensity at £ along a path direction determined by k 
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I 


I^(t; p,,cp) 

I. 

1 

L(k') 

I. (i; n,cp) 
I 

s 

Is(k) 

Ig 

sa 
I . 

SI 

I^(C',k) 

(t , k) 

Vl 

^pU 

'tL 

'rU 

I± 

J(C',k) 

^ ^33 

/< 

k 


downwelling diffuse intensity 
incident intensity 

downwelling intensity at the upper surface of a canopy 

incident intensity on a layer at optical depth t 

scattered intensity 

scattered intensity from da 

scattered portion of the intensity field at x 

intensity outgoing from 

incident intensity scattered into 

total intensity field 

total intensity field at t 

upwelling diffuse intensity 

reduced incident intensity 

radiance reflected from a lower surface 

radiance reflected from an upper surface 

radiance transmitted from a lower surface 

radiance transmitted from an upper surface 

intensity directed upward or downward 

intensity scattered (per unit length per unit density) into 
the beam at £ in a direction given by k 

unit vector; a function of 0,cp 
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A. 

k' 

A 

k. 

J 

koOo.cPo) 

A 

k 

L 

L' 

L^„(d,z) 

L^„(-d.z) 

L± 

2 

M = 1,2, 3, 4 

N 

n 

A 

n 

n da 
n da 
nJ 

nJ(£,k) 

n-kp = cos Yq 

n AV 

-nal 

P(k.k’) 


unit vector in the direction of incident intensity 
unit vector with polar angles 0^. and cp^. 

unit vector in the direction of solar or specular intensity 

A 

upwelling or downwelling k 

beam attenuation function 
lower leaf side of interest 
lower hemisphere 

uniformly diffuse (Lambertian) upwelling intensity at level 
z in the a layer 

uniformly diffuse (Lambertian) downwelling intensity at 
level z in the a layer 

uniformly diffuse upwelling and downwelling radiance 
path 

integration limits as defined in table IV 
number of elements in AV or in dV 
number density of scattering volume elements 
unit vector perpendicular (normal) to an area 
differential leaf area index 
vector area of magnitude da 
intensity scattered into the beam 

intensity scattered per unit length into the beam at £ in a 
direction given by k 

number of elements in a volume element 

intensity removed from the beam by scattering or absorption 

phase function relating incident to scattered intensity for a 
volume element of material 
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p 

Pa 

Pi 

Ps 

P± 

P± 

<P+> 

Q 

R 

R' 

S 

S(k,k') 

®L’®U 

sie^,n) 

u 

U’ 

x,y ,z 
z 

^soil 


power 

power outgoing from da into A^2 

s 

total power incident on volume element AV 

total power scattered from the volume element 

upwelling or downwelling power due to downwelling 
incident intensity 

upwelling or downwelling power due to upwelling incident 
intensity 

average power 

number determined from boundary conditions in equations 
(95) to (98) 

number determined from boundary conditions in equations 
(95) to (98) 

average value of 

surface of volume element AV 

bidirectional reflection function for the canopy as a whole 

lower and upper bidirectional reflection functions for a 
layer x' thick 

leaf zenith angle density function 

lower and upper bidirectional transmission functions for a 
layer x' thick 

positive and negative regions of ^ (t) 

upper leaf side of interest 
upper hemisphere 

coordinate frames related by rotation matrices 
horizontal canopy level 

distance from the top of the canopy to the soil 
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a 


fraction of intercepted incident power scattered 


= ±g 


P 


fraction of total power scattered into Afl^ 


r = sin 0 sin cp 

Y angle between unit vectors 

AA. cross section shared by AV and the incident beam 

AAg cross section shared by AV and the scattered intensity con- 

tained in A£2 

s 

Ap^ power incident on the volume element from intensity contained 

in the solid angle Af2^ 

Apg total power scattered into the solid angle Af^^ 


AV 

A£2 

An. 

1 

An 

s 


An 


0 



8 = cpi - 9^0 

6 .. 

1 ] 

8(p - up ,6(9 - 9+) 


volume element 
element of solid angle 

element of solid angle containing incident intensity 
element of solid angle containing scattered intensity 
solid angle that is intercepted by the Sun 

amount of Ap^ actually intercepted 

Kronecker delta function 
Dirac delta functions 


s ,e ,8 unit vectors defined in equation (120) 

X y z 

0 polar zenith angle 

0^ polar zenith angle of the unit normal n 
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V,T) 


P(k,k») 

a 

T 

A A 

T(k.k’) 


9 


9 


n 


^nO 

'I' 


cosine of the polar zenith angle in a given coordinate frame 

combined parameters indicating the functional dependence of 
the density function 

bidirectional reflection function 

soil bidirectional reflection function 

average scattering cross section of scattering centers 

optical depth 

specific optical depth 

bidirectional transmission function 

upper surface of the canopy 

polar azimuth angle 

polar azimuth angle of the unit normal n 
azimuth angle of n 

integration variable used in equation (183) 


Subscripts: 

d 

i 

i 

i 

3 

j 

L 

m 

r 

s 


downwelling 
incident direction 
summing index for p ,t 
summing index for 1,2,3 
summing index for U ,L 
summing index for 1,2,3 
lower surface incidence 
multilayer constituent index 
radial 

scattered direction 
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T 

t 

U 

u 

a 

X 

0 

1 

2 


integration interval 
total 

upper surface incidence 
upwelling 
layer index 
wavelength 

quantities associated with the solar beam 
class that is dependent on cp^ 

class that is not dependent on cp^ 

upwelling or downwelling 


Mathematical notations: 

* leaf coordinate frame 


{} 

[] 


choose the greater of two quantities contained therein 


choose the lesser of two quantities contained therein 


GENERAL APPROACH 


All analytic canopy modeling begins implicitly or explicitly with the radiation 
transfer equation that describes changes in monochromatic intensity or radiance . 
The radiation transfer equation is written 


= -n(C)d(C,k)I(8 ,k) + n(C)J(8,k) 
ax 


( 1 ) 


The equation is a statement of conservation of energy by which radiation transfer 
through any medium is governed. The change in intensity l(fi ,k) along a path ele- 
ment dfi in the direction of the unit vector k is the difference between the inten- 
sity absorbed and scattered out of the beam and the intensity scattered into the beam . 
The symbol n is the number density of scattering centers , and d is the average 
scattering cross section of scattering centers . Intensity per unit volume scattered 
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into the beam at £ in a direction given by k is represented by the expression 
nJ (£ ,k) . Integrating factors gives a formal solution 


I(£,k) 




J(C ,k)e 


-T(C ,C) 


a 


d£' 


( 2 ) 


where 


X 

no d£ 


0 


(3) 


and the J term is given by 


1 27t 


J(£-,k) - ^ If I^(£' ,k)P(k,k')dp' dcp' 


-1 0 


(4) 


The phase function P(k,lj') relates incident to scattered intensity for a volume 
element of material . The expression dp' dcp' is -dfl' where difl' is a differential 
solid angle, and p' = cos 0' where 0' is the polar zenith angle and cp' is the polar 
azimuth angle. Because the expression I^(C',k) represents the total intensity 

field at d£ , I^(C ,k) couples the solution along one path with the solution along 

all other paths; therefore, the formal solution for a single path is not the solution 
for the medium as a whole , and the solution for the medium as a whole involves an 
infinite set of coupled integrodifferential equations for which a completely general 
solution is not known . 

The problem of radiation transfer in any medium has two main parts . The 
physical problem is to find the phase function , which will be expressed in terms 
of the scattering and absorptive behavior of the individual scattering components 
or scattering centers in a volume element . The mathematical problem is to solve 
the system of coupled integrodifferential equations . 

Specification of the phase function for a given volume element will be the same 
regardless of the simplifying assumptions subsequently made to solve the mathemat- 
ical problem . That is , the phase-function development is the same regardless of 
whether parallel homogeneous layers are later assumed or whether canopy row 
effects are considered . The scattering and absorptive behavior of a volume element 
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of scattering components is a type of sum of the scattering behavior of each compo- 
nent. The problem is to assemble this sum in terms of the scattering and absorption 
properties of each component . 

The first assumption is that any individual scattering component , such as a 
leaf, may be represented by a collection of small planar elements. This assumption 
is based on the fact that any surface with any curvature can be represented to any 
desired degree of accuracy by a collection of planar elements . These planar ele- 
ments will be designated by the vector surface element da. (A vector surface 
element is a vector of magnitude da proportional to the area of da in a direction 
perpendicular outward to da.) This equation will be written 


da = dan 


( 5 ) 


where n is a unit vector normal from da in the outward direction . Further , the 
number volume density of surface elements is designated by n . For a canopy inte- 
grated over the height of the canopy , n da will give the leaf area index and thus 
can be called the differential leaf area index . 

Each element can be characterized by a slope and an azimuth in a z-vertical 
coordinate system . Because the slope i^ the same as the zenith angle of da , two 
angle parameters are associated with da: the polar zenith angle of the normal 0^ 

and the polar azimuth angle of the normal cp^ . The collection of elements in a vol- 
ume element will have some distribution of slopes and azimuths that can be described 
by density functions. The quantities volume number density, element area, and 
slope and azimuth density functions constitute the set of input parameters henceforth 
referred to as phase-function geometry . (The product of density and leaf area is 
the one-sided leaf area per unit volume.) When combined with the spatial variation 
of the phase function over the canopy, these quantities constitute the canopy geomet- 
ric inputs. These geometric inputs are a set of morphological descriptors determin- 
ing canopy reflectiAdty; they can be labeled spectral morphological descriptors. 

A planar scatterer can be parameterized by bidirectional reflection and trans- 
mission functions. For a leaf, this procedure will generally require a bidirectional 
reflection function defined on the upper surface , a bidirectional reflection function 
defined on the lower surface, and a bidirectional transmission function. These three 
bidirectional functions completely specify the scattering and absorptive behavior 
of the canopy component and are referred to in this paper as the component optical 
properties . Thus , the physical inputs to the phase function , and hence the model , 
are these optical and geometric quantities . 

The output of the model is the bidirectional scattering behavior of the canopy 
as a whole; it is expressed in terms of canopy bidirectional reflection and transmis- 
sion functions . These functions can be labeled the canopy optical properties . Once 
the canopy optical properties are known , the upwelling intensity is obtainable for 
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any specified illumination conditions from a relationship of the form 


TX 

2 27t 

I(k) = ^ // 1. (k')S(k,k')|J.' dp' d9' (6) 

0 0 


The expression I(k) represents the upwelling intensity at the upper surface of 
the canopy; I. (k') is the downwelling intensity at the upper surface of the canopy; 

/y A ^ 

and S(k,k') is the canopy bidirectional reflection (upwelling, scattering) function 
for the canopy as a whole. Thus, once S(k,k') is known for a canopy, the down- 
welling intensity field must be known to obtain the upwelling intensity . 

In summary , for operational use of the model , the input or independent vari- 
ables are the canopy geometry , the component optical properties , and the illumina- 
tion conditions . The output or dependent variable is the upwelling intensity . 

No solution is known for a medium for which the phase function varies arbi- 
trarily as a function of position . Such a variation might be a horizontal variation , 
such as row effects , or a vertical variation in which density is a function of height . 
To solve this problem , simplifying assumptions must be made . These assumptions 
may be made in terms of (1) canopy geometry, (2) component optical properties, or 
(3) first-order approximations of the field behavior within the canopy . The plant 
canopy reflection models developed in this paper vary in form and approach depend- 
ing on the combination of simplifying assumptions taken . The assumption that the 
canopy can be approximated by homogeneous layers is made in all models . Although 
the explicit development is for single-layer, single-scattering-component canopies, 
a method for generalizing to multilayer , multiconstituent canopies is given . 

Model I and Model II are closely related . The mathematical operations in both 
models are made tractable by assuming that the field within the canopy can be initially 
approximated by a uniformly diffuse downwelling component , a uniformly diffuse 
upwelling component, and an attenuated incident or specular component. The atten- 
uated specular component is defined as the portion of incident radiation at any level 
within the canopy that has not intercepted a canopy component . Under this assump- 
tion , the radiation transfer equation reduces to easily solvable modified Kubelka- 
Munk or Shuster equations . Model I is explicitly derived for leaves that need not 
have constant bidirectional scattering functions; that is, it is not restricted to the 
special situation of Lambertian leaves . However , Model II is derived solely for uni- 
formly diffuse or Lambertian leaves . Although the development techniques have 
some variations , Model II is in essence an application of Model I to Lambertian leaves 
In Model III , the most general model , the assumption is only that the canopy is rep- 
resentable by horizontally homogeneous layers . This treatment generalizes a method 
of Chandrasekhar (ref. 1) to generate a set of eight integrodifferential equations 
containing only the phase functions and the canopy reflection and transmission func- 
tions . The eight equations are reduced to six by transmission reciprocity . Model 
III terminates with these six equations , for which no general solution is known . 


13 



All three models are alike in that they begin from the radiation transfer equa- 
tion. Model I and Model II are different from Model III because a different mathemati- 
cal approach is used to solve the radiation transfer problem for the canopy as a 
whole . Identical accuracy should not be expected for all models . Model III should 
be the most accurate, and Model I should be more accurate than Model II. In partic- 
ular, Model II becomes inaccurate at wavelengths of less than approximately 0.55 
micrometer , where the assumption of constant bidirectional component scattering 
functions (Lambertian leaves) is not good. This breakdown in the high absorption 
regions is also noted by Smith (ref. 2) in the i model in which he uses this assumption. 
(Note: In a region of high absorption, scattering tends to become more specular 
because multiple scattered contributions diminish and leave primarily the more 
specular Fresnel contributions as reported by Breece and Holmes (ref. 3) .) 


FORMULATION OF THE PROBLEM 


The canopy radiation transfer problem is formulated by a discussion of the 
phase function and the formal solution to the radiation transfer equation . 


The Phase Function 

Consider a beam of cross section AA within a canopy (or any scattering 
medium) along a path £ . An element of length dC defines a cylindrical volume ele- 
ment AV. The change in intensity I in d£ along £ gives the equation of transfer 
with all functional dependence suppressed. 


d£ 


-nal + nJ 


(7) 


The term -ndl represents the intensity removed from the beam by scattering or 
absorption . The term nJ represents the intensity scattered into the beam from 
radiation incident on the volume element from directions other than along £ . This 
term could also include emission into the beam; however, emission is ignored in 
this discussion. 

The quantity n is the number density of scattering components in AV of 
a given type . Then , the fraction of the beam scattered or absorbed will be given 
by the fraction of AA obscured by the components in AV . The number of the com- 
ponents in AV is n AV, and k is the beam*attenuation function. Thus, the 
fraction intercepted is 


, n AVa 

‘ Ta- 


( 8 ) 
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where a is the average cross section of components in the direction E . However , 
along a cylindrical path 


AV 

AA 


= d£ 


(9) 


Thus , for the change in intensity due to intensity removed from the beam 

^ = -nol ( 10 ) 

A particular form for the average cross section a will be developed subsequently . 
Meanwhile , it is sufficient to note that a will be a function of the direction of £ , 
a complication not encountered in the simpler problem of atmospheric radiation 
transfer . 

An expression for the term J will also give an expression for the very impor- 
tant, phase function . The expression nJ denotes the radiation scattered into the 
beam. To derive this expression for nJ , consider an arbitrary volume element AV 
having a cross section AA shared by both the beam and the scattering medium . 

The total power incident on this volume element p^ is found by integrating the de- 
fining equation for intensity over S , the surface of AV . This integration is written 


4n 

Pi = f y'L(k')(k'*di)df2. 

0 S 


( 11 ) 


A 

where k' is a unit vector in the direction of incident intensity I. and where ds 

A ^ 

is a differential vector area on the surface of AV. Note that k is a function of 9 
and 9 ; therefore, a function of k and k' is also a function of 0,cp;9',cp'. The 
power incident on AV from the intensity contained in Afl^ , an element of solid 

angle containing the incident intensity , is designated Ap^ and is written 


Apj = L Af2j^^k'*ds 


( 12 ) 


But 


fi' 


ds = AA. 

1 


(13) 
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where AA^ 


is the cross section shared by AV and the incident beam. Then 


Ap. = I. An. AA. 
^1111 


( 14 ) 


Similarly , the total power scattered into the solid angle An^ is denoted by 

Ap and is written 
^s 


Ap 


s 


= I AH AA 
s s s 


(15) 


where AA is the cross section shared by AV and the scattered intensity con- 

® 2 
tained in An . The amount of Ap. actually intercepted, A p. , is written 

S 1 1 


2 _ AV 

Pi = ^‘^aa: 




(16) 


If the fraction of 


a 2 

A p, 


scattered is a, then 


P 


s 


.2 

a A p. 


(17) 


Furthermore , if the fraction of p scattered into An is denoted by p , then 

s s 



An 

s 



(18) 


Using equations (14) , (16) , and (17) , the following is obtained. 


^Ps 


An 

s 


= apno 


AV 

AA. 

1 


I. An. AA. 

Ill 


(19) 
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The expression Ap is obtained from equation (15) . Thus 

s 

I . 

— = aBna ACl.l. (20) 

) 

where the subscript i on I denotes the scattered intensity resulting from the 
scattering of L contained in . Let P(k,k^) be defined by 



P(k,k^) = 47rapo 


( 21 ) 


Then 


P(k,k^) 


47tI .(k) 
SI 


n 


AV 

AA 


I.(k')An. 
1 1 


( 22 ) 


The phase function P(k,k^) relates incoming and outgoing intensity for an 
arbitrary volume element. This formulation of the phase function is very useful 
in developing the form of the phase function for the canopy . Note that , in general , 
the total scattered intensity from a volume element has contributions from incident 
intensities in all solid angles. Summing equation (20) over all AH. values yields 
an integral ^ 



471 

i.(k')P(k,k')da' 

0 


(23) 


where I (k) is the total scattered intensity in the k direction as a result of scat- 
s 

tering of incident intensity by the volume element AV. When the expression for 
a cylindrical volume element (eq. (9)) is recalled, it is clear from the discussion 
of equation (1) that equation (23) is also the equation for n J . Thus, written in full, 
showing all functional dependence , the expression for J is 


471 

J(C,k) = I^(C,k')P(k,k')da' 

0 


(24) 
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Formal Solution of the Radiation Transfer Equation 


The radiation transfer equation is 


dl (£,k) 
dfi 


= -n(e)a(e,k)I(E,k) + n(2)J(fi,k) 


( 1 ) 


T(2,2q) 

This equation can be integrated by multiplying each term by e where 

T(£ ,£q) is defined in equation (3) . Note that from a rule of definite integrals 


and 



(25) 



T(£,£') 



(26) 


where 2^ < £' < 2. Then, both sides may be integrated from £ to 2^ to obtain 
equation (2) . 


1(2, k) 


■.M.-’C'-)./ 


J(2' ,k)e 


-T(2,2') 


d£’ 


( 2 ) 


with the quantity J (2' ,k) given by equation (4) . 


The formal solution gives intensity at 2 due to scattering along the beam 
before 2^. If the value of 2 is at a surface of the scattering medium, the solution 

gives scattered intensity from the medium . For a canopy , this could be the quantity 
desired, the upwelling intensity. The problem consists of finding the right-hand 
terms. The first right-hand term is simply the intensity at some point 2^ on 2. 

This intensity might be determined by boundary conditions and , as such , may not 
be an insurmountable difficulty. The real problem arises from the second right-hand 
term and in particular from the integral for J . 
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Expressing J requires knowledge of the phase function. As will be shown, 
to obtain a canopy phase function is difficult but not impossible . However , the J 
integral also requires knowledge of the value of I at any point along 2 . If this 
value were known , there would be no need for the formal solution . Thus , there 
is an infinite set of coupled integrodifferential equations , a computational problem 
for which there is no direct solution known to the author . 

A canopy phase function that is based on a minimum of assumptions and that 
is believed to be completely general within these assumptions will be developed in 
this paper. Two techniques for handling the computational problem will be explored; 
both techniques begin from the assumption that the scattering medium can be repre- 
sented by parallel layers . For a canopy , this assumption will mask row effects . 

One of the computational methods , Model III , a modification of a technique developed 
by Chandrasekhar (ref. 1) for atmospheres, involves restatement of the problem 
in terms of other system parameters , namely reflection and transmission functions 
for the layer as a whole . Then , the radiation transfer equation and its formal solu- 
tion are used to solve the parameterized problem . This method concludes with a 
finite set of integral equations . The solution to these equations may provide ulti- 
mately the most exact technique for developing a general canopy model. In the 
sfecond computational method , used for Model I and Model II , the modified Kubelka- 
Munk equations are used to find a first approximation to the intensity field within 
the canopy . This technique is simpler and may be found workable from a practical 
computational point of view . Improvement of the Kubelka-Munk approximation can 
be obtained by iteration . 


BIDIRECTIONAL PHASE FUNCTION FOR A CANOPY 


The so-called phase function shown in equation (22) would be better named 
volume scattering- absorption function; it is basically a ratio of the incident to the 
scattered intensity for an arbitrary volume element within a medium . The scattering 
behavior of the volume element is a sum of the scattering behavior of the individual 
components , or scattering pieces , in the volume element . 

Although performing this sum is straightforward in principle , it is fairly 
complicated in practice. The general summing problem is expressible as the integral 
sum of sectionally continuous functions in four dimensions. Each dimension is non- 
zero in half spaces , defined with respect to coordinate frames , each of which is two- 
dimensionally rotated with respect to the other . The solution of this summation prob- 
lem , the coordinate transformations , and the integration limits are important at this 
stage of theory development . 


Assumptions and Canopy Visualization 

The canopy is viewed as a collection of volume elements , each of which con- 
tains an ensemble of canopy components. Each component, such as a leaf, is as- 
sumed representable by one or more planar elements. Each planar element, repre- 
sented by the vector area da, is characterized by slope, horizontal orientation or 
azimuth, and optical scattering properties described by upper and lower surface bi- 
directional reflection functions and a bidirectional transmission function . The slope 
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is the zenith angle in a z-vertical coordinate frame of the unit normal to da, and the 
azimuth is the azimuth in this frame of the unit normal to da . The zenith and azimuth 
distribution for the ensemble of elements in the volume element will be characterized 
by statistical density functions . These angular density functions , together with the 
volume number density and the size of area elements, will be called the plant geom- 
etry. The bidirectional reflection and transmission functions will be called the com- 
ponent optical properties , The plant geometry and component optical properties 
are the necessary inputs to the model. 

Complete optical properties of canopy components are available for only a very 
limited number of components . These properties have been measured satisfactorily 
only by Breece and Holmes (ref. 3) for corn and soybeans. Therefore, in current 
practical usage , estimates must be made from the numerous directional reflection 
measurements. This procedure is equivalent to assuming that the components are 
uniformly diffuse or Lambertian . Analysis of available data suggests that such an 
approximation is acceptable for wavelengths greater than approximately 0.7 microm- 
eter. However, for wavelengths less than 0.7 micrometer, the model will not be 
expected to work well without a better estimate of the nonuniform reflection proper- 
ties of the leaves . This result is borne out by the modeling work done by Smith 
(ref. 2) . 

Polarization will be disregarded in this model . However , polarization is an 
important effect that will likely prove to be a useful diagnostic tool . To account for 
polarization , the formal theory presented herein can be used with the understanding 
that the scattering functions are thereby represented in polarization matrices . The 
aspect of polarization most significant to remote sensing is the existence of a strong 
front surface quasi-specular reflection for many leaves , notably waxy ones . This 
reflection or gloss is not greatly affected by the absorptive character of the leaf. 

For many aspects of spectral signature interpretation, this gloss masks the spectral 
information and could be called noise . The same effect is seen in an ordinary glossy 
photograph viewed in glaring light. This glare is often strongly polarized; thus, 
polarization filtering could be used to remove it. For this reason, polarization re- 
search is potentially valuable and could make an important contribution to remote 
sensing . 


Coordinate Frames 

The development of the canopy function requires the definition of three coor- 
dinate frames and the transformations between them. In the first frame, the obser- 
vation frame, the z-axis is vertical, the xy-plane is horizontal, and the x-axis is 
left arbitrary to be chosen in any way that may become convenient later . In the 
second frame, the upper leaf frame denoted by asterisks, the z*-axis is perpendicular 
to the planar element; thus, the z*-axis has the same direction as the unit normal 
n to the element da. The x*-axis is chosen along some convenient line of the com- 
ponent, such as the central vein of a leaf. Then, the x*y*-plane is the plane of the 
element da. The third frame is the lower leaf frame denoted by double asterisks. 
The upper and lower leaf frames are related by z* = -z**, y* = -y**, and x* = 

X** . Thus , the double-asterisk frame is related to the asterisk frame only by the 
direction of the outward-drawn normal . 
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Each of these frames is conveniently expressed in spherical polar coordinates 



using the convention where 0 is the zenith angle and cp is the azimuth angle meas- 
ured counterclockwise from the x-axis when viewed from the upper plus-z hemi- 
sphere in the minus-z direction. The assumption is that the planar element is char- 
acterized by some slop^nd azimuth. These are the zenith and azimuth angles of 
n, the unit normal to da, measured outward from the upper surface. Let 0 ^ be 

the zenith angle and be the azimuth angle of n in the observation frame. 

The transformation equations are obtained by first considering a rotation of 
9^ around z by which y is moved to y* , and then taking a rotation of 0^ 

around y*. The sum of these two rotations will take z into z* and x into x*. 

The first rotation matrix around z is 


C = 


cos 9^ sin 9^ 0 


-sin 9 cos 9 0 

^n ^n 


( 27 ) 


The second rotation matrix around y* is 


V = 


cos 0 0 -sin 0 

n 


n 


sin 0 0 cos 0 

n 


n , 


( 28 ) 


The product is 


R = VC = 


cos 0 cos 9 cos 0 sin 9 -sin 0 


n ' n 
-sin 9^ 


n 

cos 9 


n 


sin 0 cos 9 sin 0 sin 9 cos 0 


n 


n 


n 


n 


n 


( 29 ) 
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Now, consider some radial unit vector in some direction k^. Expressed in the ob- 
servation frame , this term may be written 




sin 9 cos q) 
sin 0 sin cp 
cos 6 


\ 

/ 


The same vector expressed in the asterisk frame is 


k* 

r 




sin 0* cos cp 
sin 0* sin cp 
cos 0* 



But 



(30) 


(31) 


(32) 


Applying the rotation matrix to k^ , comparing components with k* , and solving 
gives the transformation equations 


-cp** = cp* = cot 


-1 


cos 0 


n (‘P ■ ‘Pn) 


sin 0 cot 0 
n 


sin ^cp - cp^^ 


(33) 


cos (7t - 0**) = cos 0* = 


cos 0 cos 0 + sin 0 sin 0 cos 
n n 


(q> - 


-p.** = p* = COS 0^ cos 0 + sin 0^ sin 0 cos ^cp - 9^^ 


(34a) 


(34b) 
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Bidirectional Reflection and Transmission Functions 

The derivation of a phase function for an ensemble of canopy components con- 
sists of two main parts , The ^st part is to relate incident to scattered intensity 
for an individual component da , and the second part is to sum the contributions 
of all the components in the volume element . The reflection and transmission bi- 
directional scattering functions are defined by equations of the form 

1 2n 

Ig(k) = I y'y* p(k,k')I.(k')ji' dp' dcp' (35a) 

D 0 

where 0 < p < 1 , and 

1 2n 

Ig(k) \ f f t(k,k')l.(k')p' dp’ d9' (35b) 

0 0 

A. A. A 

where 0>p>-l and where p' arises from k'»n wiftin n , a unit normal to da. 
The expression 1 (k) is the scattered intensity from da, and I.(k') is the incident 

® A A A ^ 

intensity on da. The quantities p(k,k') and t(k,k') are the bidirectional reflec- 
tion and treinsmission functions , respectively . An alternate form could also be used . 

1 2iz 

Igfk) ^ ^ f f P(k.k')l.(k')dp’ dcp' (36a) 

0 0 

where 0 < p < 1 . 

1 2x 

*s^^^ "" 4^ y * f ^(^’k')Ii(k')dp' dcp' (36b) 

0 0 

^ A 

where -1 < p < 0 and where p = k«n. Note that equations (36a) and (36b) are 
obtained from equations (35a) and (35b) by omitting the p' within the integral and 
inserting the l/4p before the integral. In both equations , the bidirectional func- 
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tions so defined are symmetric under interchange of directions . That is 


p(k,k’) = p(-k’.-k) 


(37a) 


and 


T(k,k’) = T(-k',-k) 


(37b) 


Because da is a planar element that will generally have different reflective 
properties depending on the side of incidence , it is convenient to define the total 
reflection and transmission functions as shown in table I. All the bidirectional func- 
tions shown in table I are zero in regions other than those indicated. Note that limits 
are given in the asterisk or upper leaf frame where the plane of the leaf is the xy- 
plane. The subscripts U and L designate upper and lower sides of incidence, 
respectively . Radiation scattered into the hemisphere of incidence and radiation 
scattered into the hemisphere opposite the hemisphere of incidence are governed 
by the transmission function t(k,k’). 

In general , the scattered field will be given by 


Ig(k) = Ipu(k*) + (38a) 

where 0 < p,* < 1 , and 

Is(k) = 


where -1 < u* < 0 and where I 


I 


and I 


are derived from the equa- 

The 

existence of a transmission-function symmetry relationship that makes the distinc- 
tion between and redundant can be shown . However, the upper and lower 


P Xpy, 

tions defining the bidirectional scattering functions (eqs. (35a) and (35b)) . 


surface transmission-function symmetry will not be explicitly exhibited in the deri- 
vation because omitting it makes checking the derivation less complex . Explicitly , 
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this derivation is written 


1 2n 

Ipy(k) = ^ // I.(kOlk'*nlPuCk*,k’*)d^' d9» (39a) 

0 0 

1 27t 

I^U<k) = ^ // L(k')|k'»n| Ty (k*,k’*)d|a' dcp' (39b) 

0 0 

0 2n 

// I. (k') |k '‘nl Pj^(k*,k**)d^i' d 9 ' (39c) 

-1 0 




1 

Tt 


0 27T 

// I.(k') |k'*n|Tj^(k*,k'*)dia' dcp’ 
-1 0 


(39d) 


Note that no terms in these expressions have coordinates with asterisks ex- 
cept the bidirectional scattering functions . The necessity for expressing all scat- 
tered intensity in the observation frame (z-vertical) ultimately, even though the 
bidirectional scattering functions are customarily measured in the leaf frames (with 
asterisks) , is emphasized by this fact. The absolute value on the dot or scalar prod- 
uct in the integral ensures only positive values of reflected intensity . Note that 
the vector intensity notation indicates the specific p. included in the integral . Re- 
member that intensity may be conveniently expressed as a vector for some purposes; 
however , it does not always follow the laws of vector addition nor is it a vector field 
in the same sense that an electrostatic field is a vector field . The intensity field has 
a magnitude for every direction associated with every point . The electrostatic field 
has only one magnitude and direction associated with every point. Further, two 
beams do not add vectorally for purposes of reflection . 

Because any incident radiation field on da can be viewed as a collection of 
beams contained within nonoverlapping , small solid angles , the basic problem can 
be reduced to considering incident radiation in some solid angle . Integration 
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gives the two pairs 




if and only if incidence is on the upper surface , and 


(40a) 


(40b) 


(40c) 


(40d) 


if and only if incidence is on the lower surface, where the nonzero regions are de- 
fined in equations (38a) and (38b) . The dot product of a unit vector k in a direc- 
tion of interest j and a unit vector n normal to da expressed in the observation 
frame without asterisks gives 


k.*n = cos Y- = sin 6. sin 0 cos /cp. - cp \ + cos 0. cos 0 (41) 

] ^3 ] n 1^3 ^nl 3 n 


where y is the angle between unit vectors . By applying the transformation equa- 
tions and integrating, the complete scattered intensity field for the individual ele- 
ment da is obtained. Note, however, that this field is dependent on the polar angles 

0 and cp . 
n ^n 


Canopy Geometric Density Functions 

The next step is to find the intensity of an ensemble of planar elements , each 
of which is characterized by a slope and an azimuth. Let and Ng describe the 

number of elements da in the volume element AV that have a normal between 

cp and cp + Acp and 0 and 0 + A0 , respectively . Let N equal the total 

^n ^n ^n n n n ^ 
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number of elements in AV , Then , density functions may be defined by 



lim ^^^n 

Acp ->0 N Acp 
^n ^n 


lim 

A0 -»-0 
n 


AN0 

n 

N A0 

n 


( 42 ) 


(43) 


where v and r) are combined parameters indicating that D(cp^,v) and S(0 ^,t]) 

may depend on such factors as stress , age , wind , and diurnal effects . The density 
functions D(cp^,v) and S(0^,r|) then compose the fraction of leaves in the inter- 
vals cp to CD + Acp and 0 to 0 + A0 and thus are subject to the normali- 

^n ^n ^n n n n •’ 

zation condition 


27t 


/ 


D 



1 



(44) 


(45) 


The upper limit on S(0^,r|) implies that the normal to the upper surface of da 

never drops below the horizontal . This simplifying assumptjnn is convenient but 
not necess^ary. If this assumption is not true, then values of da having a downward- 
pointing n constitute another class of da that must be handled similarly to those 
not having a downward-pointing n . 


Formal Summing of Planar Ensemble Scattering 

Having calculated the scattered radiance ^r an individual element da , the 
problem is to sum the scattered radiance of all da in volume AV. Equations (38a) 
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and (38b) can be compactly written 


I 

sa 


= ? 'ija 
IJ ^ 


(46) 


where i = p,t and j = U,L. The defining equation for intensity is 


1 2n 

p = // I(k)|k*daldag (47) 

0 0 

Therefore , the power outgoing from da into Afi is 

s 


where I (k) is intensity outgoing from da into . The total power into AQ, 

S& S 1 

from all da in AV is 


Ps = Pa 

a=l 

This sum over all p is accomplished by using the density functions D(<p ,v) 

8. n 

and S(0 ,n). The power into AQ of all da in AV in intervals d0 and dcp 
n ' s n ^n 

at 0^ and cp^ then is given by 


dp^ = n AVD^cp^.vWe_^,t,^ da-I^^ A£2^ d<p^ de^ 


Thus 


Ps ' n AV Afi^yi^^(k)|k^.dalD^(p^,v^S^0^,ti)dcp^ d6_^ 


(50) 


(51) 
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t 


where the integration limits over and 0^ must be carefully defined . These 
limits will be considered later . 

From equation (15) , the power scattered from AA into Afl may be ex- 
pressed by 


P - I 

^ c* < 


{K) 


AH AA 
s 


(52) 


Equating expressions for p gives 

s 




(53) 


Using equation (22) , which defines the phase function, P(k,k') can be written 


P 



“’"s (\) 
AA i\ ij 1 


(54) 


Then 




_ 4an/p(<p^.v)s(e^,n)d9_^ 


1. Af2. 
1 1 


(55) 


However, from equations (38a) and (38b) 


I = I ,, + 1 ,, + 1 ^ + I , 
sa pU tU pL iL 


(56) 


Note that 


I 

sa 



= I 

sa 



da cos 


(57) 
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Then, using equations ( 39 a) to ( 39 d) and ignoring integration limits, P(k ,k.) 
may be compactly written ® ^ 


P 



/ 


47m / COS Y* cos Y ip 




u 



X D 



d0 


n 


( 58 ) 


where cosy- is a function of 0.,cp.;0 ,cp , and cosy is a function of 0 ,<p ; 

1 11 Ti- n s s s 

0 ,cp , Transformation to the observation frame makes the bidirectional functions 
n n ^ ^ 

dependent on ®g> 9 g> Note also that P(k,k') is symmetric under 

interchange of directions . 

The phase function has been essentially reduced to an integral , the limits of 
which are not straightforward or obvious and must be carefully defined . The indi- 
vidual bidirectional reflection scattering functions for an_ individual element are 
noted as potentially different depending on the side of da on which the radiation 
is incident. In addition, the transformation of the bidirectional scattering functions 
from the observation frame with asterisks to the observation frame without asterisks 
has a complication . These functions are defined only in half spaces separated by 
the plane of da. Therefore, the integration limits must be carefully specified to 
cover only the nonzero regions and thereby to prevent the appearance of contribu- 
tions from noncontributing regions . 


Phase-Function integration Limits 

For a given incident and scattered direction , the problem is to specify the 
integration limits on equation ( 58 ) so that integration includes only domains for 
which the particular bidirectional function is nonzero in the z-vertical observation 
frame. The quantity Ip^ will be nonzero for leaves both illuminated and viewed 

from the upper surface; Ipj^ will be nonzero for leaves illuminated and viewed 

from the lower surface . The quantity 1 ^^ will be nonzero for leaves illuminated 

on the upper surface and viewed from the lower; will be nonzero for leaves 

illuminated on the lower surface and viewed from the upper. The nonzero radiance 
for a given leaf depends on the combination of 0 and 9 ranges that character- 
izes the leaf. ” ” 

The fundamental calculation used to define these various domains is calcula- 
tion of the conditions for which a direction of interest is incident on the upper sur- 
face or incident on the lower surface . This direction of interest will be either the 
illumination or the view direction . This calculation may be done by considering 
Wie dot or scalar product between the outward normal from the upper surface of 
da and a u^t vector in^ the direction of interest. Let n be the unit vector perpen- 
dicular to da and let k. be a unit vector in the direction of interest. The dot 
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product between unit vectors gives the cosine of the angle y between them. This 
quantity, which was shown in equation (41) for y^ , is 


cos y. = sin 0. sin 0 cos ( cp. - cp \ + cos 0. cos 0 (59) 

'i in ^n) 1 n 


If cos y. is positive, incidence is on the upper surface; if it is negative, in- 

cidence is on the lower surface . When cos y^ is zero , n is perpendicular to 

and, thus, the planar element is parallel to the direction of interest. When 
cos y. = 0 


cos 6 = -cot 0. cot 0 
1 n 


( 60 ) 


with 



(61) 


where cp^^ is the azimuth angle of n at which cos y^ changes signs . If 
|cot 0. cot 0^1 > 1, then 6 cannot exist and cos y^^ does not change signs. Be- 
cause tan 0^ is considered always positive , this condition may be written 


cot 0.1 > tan 0 
1 ' n 


(62) 


which holds if 


0 < 0 

n 



(63a) 


or if 


0 < 0 

n 


< 0 . 

1 


ZL 

2 


(63b) 
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Under the conditions of equation ( 63 a) , cos y. is positive; of equation ( 63 b) , nega- 

A 

tive. Thus, in the regions defined, k. is incident on the upper and lower surfaces, 
respectively, independently of 9^. 

Requiring that cos y. be positive gives 

cos ('Pi - %) > cos 6. ( 64 a) 

where 0 < 0 . < 7t/2 or 
1 

- 9 \ > -cos 6. ( 64 b) 

i nl 1 


cos 19 




where tc /2 < 6^ < tc. Under these conditions, cos Yj is positive if 

CD. + 6 . < o) < CD. + 27 t - 6 . ( 65 a) 

* 

where 0 < 0 . < 7t/2 and 

9. + 7T - 6. < 9 < 9. + 71 + 6. ( 65 b) 

1 ^n ^11 


where 7 i/ 2 < 0 . <7t. 


Similarly , requiring that cos y . be negative gives 


('Pi ■ 'Pn) 


COS f9. - 9_1 < cos 6. 


(66a) 


where 0 < 0. < 7t/2 and 


cos 


(•Pl - 'Pn) 


< -COS 6. 

1 


(66b) 
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r 


where %/!<&.<%. Thus 
1 

cp. + 6. < cp < cp. + 271 - 8. (65a) 

^1 1 ^11 


where 0 < 0 . < ti /2 and 

cp. + 71 - 8 . < cp < cp. + 71 + 8 . (65b) 

1 ^1 1 


where 7 t /2 < 0^ < 7 t. The value of 8 . ranges from %/2 to 7 t. 

Equations (63a), (63b), (65a), and (65b) enable construction of table II , 
which gives six combined ranges of 0 ^ and cp^ for incident beams and a similar 

number for scattered beams. These combined ranges are organized by subsets de- 
pending on whether the direction of interest is in the upper or lower hemisphere 
with respect to the horizontal . The primed classes indicate lower hemisphere direc- 
tions of interest. The symbols U and L indicate whether the incident or scattered 
beam lies in the upper or lower hemisphere with respect to the plane of the leaf. 

The subscripts i and s distinguish incident from scattered beam directions . The 
subscripts 1 and 2 distinguish ranges that are unrestricted in 9 ^ from those 

that are restricted in 9 . Thus , class U „ would give the range of 0 and 9 

n sz n n 

for which the scattered intensity emanates from the upper surface of the leaf and 

lies in the upper hemisphere with respect to the horizontal. 

For a given incident and scattered direction , a leaf can lie only in one inci- 
dent class and in one viewing or scattered class . The intersection of the six illumina- 
tion classes with the six viewing or scattered classes then provides the integration 
limits over the bidirectional functions for all possible combinations of illumination 
and viewing. These intersections are labeled as classes 1 to 36 in table III, in which 

the braces | } mean that the greater of the two quantities contained within should 

be taken , and the brackets [] mean that the lesser of the two quantities contained 
within should be taken. 


Canopy Phase Function 

The phase function , including integration limits , can now be written using 
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equation (58) and table III. 



( 67 ) 


where the limits a , b , c , and d for the expression M are given in table IV . 

The entries refer to the numbers of the intersection regions in table III . The factors 

cos y. and cos y are obtained from equation (41) and may be written explicitly 
1 s 


cos y. = sin 9. sin 0 cos 
* 1 1 n 



+ cos 0. cos 0 
1 n 


(59) 


cos 


sin 0 sin 0 cos 
s n 



<p 1 + cos 0 cos 0 
^nl s n 


( 68 ) 


Note that i and s refer to incident and scattered directions for a given volume 
element , not for the canopy as a whole . These equations complete the formal speci- 
fication of the phase function. 


MODEL I - SOLUTION OF THE RADIATION TRANSFER EQUATION USING THE 

KUBELKA-MUNK APPROXIMATION 


One approach to solving the radiation transfer equation begins with a tech- 
nique for approximating the field within the canopy . When approximated , this field 
can be used to solve the radiation transfer equation using the full anisotropic canopy 
phase function . The assumption that the field can be factored to first approximation 
into an upwelling uniformly diffuse component, a downwelling uniformly diffuse com- 
ponent, and a reduced incident flux (remaining unscattered or absorbed incident 
flux) is made in this approximation technique . This assumption is equivalent to 
initially assuming a phase function that scatters uniformly across the upper hemi- 
sphere and uniformly across the lower hemisphere. The upwelling, downwelling, 
and reduced incident fluxes can be unambiguously expressed in terms of irradiance 
(power per unit area) on a horizontal plane at some level z in some layer i . Thus , 
Ej^^(d,z) is read as irradiance at z in layer a resulting from the uniformly dif- 
fuse upwelling intensity. The expression Ep^^(-d,z) is read as irradiance at z 
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in layer a resulting from the downwelling uniformly diffuse intensity. The expres- 
sion E. (s,z) is read as irradiance at z in layer a resulting from the reduced 

A(X 

incident specular flux. Each irradiance depends on wavelength X. 

Using L. (d,z) (L as in Lambertian) for uniformly diffuse upwelling inten- 
Act 

sity (radiance), L^^(-d,z) for uniformly diffuse downwelling intensity , and 

I. (s,z) for reduced incident intensity, the relationship between intensity and ir- 
Act 

radiance is given by 


E^„(d.z) = r.L^„(d,z) (69) 

E^^(-d,z) = 7tLp^^(-d,z) (70) 

E^^(s.z) = (71) 


where is the solid angle subtended by the source of the specular intensity . 

s 

No confusion should arise from subsequently suppressing the subscripts X and 
a . Note that the specular intensity could also be expressed by Dirac delta functions 
if desired . 

Under the two stream plus specular and semi-isotropic phase-function assump- 
tions , the radiation transfer equation can be transformed into modified Shuster or 
Kubelka-Munk equations . This transformation gives a set of equations 


DE 




j=l 


a..E. 
U 1 


(72) 


where i = l,2,3, D = ; 3 — , a., represents expansion coefficients, and E. and 

QZ IJ 1 

Ej refer to E(d,z), E(-d,z),and E(s,z) as appropriate. 


E [a„ - 6..D] Ej - 0 <73) 


35 



where 5„ is the Kronecker delta function . This equation has solutions if and only 
if the determinant 


det a.. - D 6.. = 0 
1 ] 1 ] 


(74) 


This equation gives the eigenvalue equation 



(75) 


Solving for gives 

k z 

F = F (0)e “ (76) 

a a 

The symbol F^ represents eigenfunctions of the D operator with eigenvalues 
The may then be expressed as linear combinations of F^ . 


3 


0=1 


A. F 
10 o 


(77) 


where expansion coefficients A. are functions of expansion coefficients 
the boundary conditions . 


a.. 


1 ] 


and 


The triplet of equations summed in equation (72) indicates that the rate of 


change in E 
cients . 


is linear with that of E.. The a., are mathematical expansion coeffi- 


1 3 

The physical task is to identify the 


1 ] 


a.. 

1] 


with physical parameters of the 


system that are measurable and that might vary with conditions to be detected. 
For example , if a^ values vary with leaf slope , which may be a function of plant 

moisture stress, then moisture stress can possibly be detected by looking at the 
ay . The physics of this approach is contained almost exclusively in the a^ . 


Formal Solution of the Kubelka-Munk Equations 

Explicitly, the Kubelka-Munk differential equations as given in equation (72) 
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can be written 





( 78 ) 


dE_ 

"d^ “ ”®21^+ ^22®- ~ ®23^s 


(79) 


dE 

~dz ^ ®33®s 


(80) 


where E = E , E_ = E_’, and E„ = E . Equation (78) indicates that in the plus-z 
direction , the change of irradiance E^ arises from E^ either scattered or absorbed 
into E_ , which is a negative change requiring the minus sign , or from E_ and 
E scattered into E . Equation (79) is similar except that the change is positive 

S ' 

in the minus-z direction; thus, the signs are changed. Examination of equation (80) 
shows that specular radiation once absorbed or scattered is no longer considered 
specular; thus, =a„„ = 0. Remember that E changes along the beam path , 

which is not necessarily along z . This characteristic will be noted in the explicit 
calculation of a^^. 

Solutions of the secular equation give the operator equations 


- ^33 = ° 


(81) 


D = a^D = a 


(82) 


where 


= 


22 


- a 


11 


* Vfcu 




^^ 12^21 


(83) 
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Using equation (77) , the following can be written. 


a^z a_z a„„z 

= A^^e + ^13® 

(84) 

a^z a_z a^gZ 

E_ = A2j^e + A22© + ^23® 

(85) 

^33^ 

(86) 

Substituting these expressions into equations (78) and (79) and equating co- 
efficients of like exponentials gives 

(“+ ^11^^11 " ®12^21 "" ® 

(87) 

^ll)^12 " ^12^^22 "" ® 

(88) 

^®33 ^1^^13 " ^12^23 " ^13^33 ” ® 

(89) 

(% " ^22)^21 ^ ^21^11 = ® 

(90) 

(®22 " “-)^22 ” ®21^12 ^ ® 

(91) 

(®33 " ^22)^23 ^21'^13 ^23'^33 " ® 

(92) 
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Solving equations (89) and (92) gives 


. _ ^21^3 ^ ^ 23(^33 ^ ^l) 

23 ^11^21 (®33 ^ ^11) (^33 ^ ^22) 

^ ^12^23 ~ ^13 (^33 ~ ^22) . 

^12^21 (^33 ■ ^22)(^33 ^ll) 


Solving the remaining four equations gives 


where 


All = (1 - <+)9 
Ai2 = (1 + <_)r 

^1 = (1 " 

^22 ~ ^ ~ 



+ a 


11 


+ a 


11 


12 


+ a 


12 




^22 ^21 
‘ ®22 “ ^21 


( 93 ) 


(94) 


(95) 


(96) 


(97) 


(98) 


(99) 


( 100 ) 
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The quantities Q and R are determined from the boundary conditions . That 
radiation is continuous across layer boundaries is assumed , If the only external 
radiance source is without a diffuse component, then 


E (0) = 0 


( 101 ) 


Further , the upwelling radiation field at the soil is given by 


E 


+ 







1 2tz 


0 0 ^ ' 


( 102 ) 


The quantity from equation (86) is clearly the value of specular radiation at 

z = 0 , the top of the canopy . The distance from the top of the canopy to the soil is 
^soil' reflectivity of the soil is Thus 


"^33 = 


<103) 


The expression a^^. = a^^ will be shown to be equivalent to assuming that the 

two sides of the scattering components (leaves etc.) are optically indistinguishable 
and are uniformly diffuse (Lambertian) scatterers. Under this assumption, the fol- 
lowing equations are obtained where is defined as ±g. 


g = 



(104) 


where a = = a^^ 


and b = a ^2 = ^21- 
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(105) 


_ c' b - c(k - a) . 

^13 2 , 2 ^33 

g - k 


where c = ~ ^23’ ^ “ ^33' 


A - cb + c' (k + a) . 
^23 2 _ , 2 ^33 




= i 



A 


11 



Q 


A 


12 



R 


Aji = (i m) Q 



The solution to the differential equations from equations (84) to (86) 
written 


E^(z) = Q(1 - + R(1 + Oe 

E_(z) = Q(1 + + R(1 - 


(106) 


(107) 


(108) 


(109) 


( 110 ) 


( 111 ) 
can be 


( 112 ) 


(113) 
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(114) 


kz 


Generalization to multilayers is accomplished by using the unabbreviated notation 
for E (eqs. (69) to (71)) and subscripting all expansion coefficients a., and A., 
and their derivative quantities with the layer label a . i] 


Development of Expansion Coefficients ^12' ^22 

The first step in deriving expansion coefficients and 

is to consider the interaction of a single planar element with a downwelling or up- 
welling uniformly diffuse field. An edge-on view of such an element is presented 
in figure 1 . The symbols and L_ designate the upwelling and downwelling 

uniformly diffuse radiance of the field . These are related to the upwelling and down- 
welling irradiance evaluated on a horizontal plane by 


irL^ 


= K 


(115) 


Let Ipy(k) be the intensity reflected from the upper surface of da; let 
1 ,, (k) be the intensity incident on the upper surface transmitted outward from the 

A 


lower surface; and let I , (k) 

pL. 


and (k) be the reflected and transmitted intensity , 


respectively, for intensity incident on the lower surface. Some portion of these re- 
flected and transmitted intensities will be upwelling and some^portion downwelling. 
These intensities are expressed as a function of the direction k, which may be ex- 
pressed by the usual 0 and 9 polar coordinate angles in any chosen coordinate 
frame . 


In this coordinate-frame-independent notation , the bidije^ctional reflection 
and transmission functions may be written Pj(k,k') and t^(k,k') with j=U,L. 

The symbols U and L refer to the upper and lower surface of incidence , respec- 
tively. Thus , the defining equations for Pj(k,k’) and (k,k') may be written 

Ip.(k) = J y* I.(k’)pj(k.k')k'*n 6Q' (116) 

M 
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and 


I .(k) = k A(k')T.(k,k’)k'*n dfl' (117) 

71 ^ 1 J 

M 

The quantity L refers to incident intensity , and integration is over the upper or 
lower hemisphere as defined by the plane of da . 

Although the subscripts on t, U and L, will be retained as a simplifying 
measure, note that reciprocity requires the transmission function obtained for upper 
surface incidence to be related to the one obtained for lower surface incidence by 
simple interchange of direction. That is, because 


T(k,k') = t(-k',-k) (37b) 


and because k 


A 

and k' 


always lie in opposite hemispheres for transmission, then 


Tj^(-k',-k) = Ty(-k’,-k) 


(118) 


Clearly, the objective of this development is to obtain the total upwelling and 
the total downwelling power resulting from the scattering of and L_ by the 

single planar element da . The coordinate frame for performing the integrations 
may be chosen in any way that simplifies the calculation . The coordinate frame used 
in this portion of the model derivation is illustrated in figure 1 and defined in table 
V . In this frame , the x-axis is down the line of maximum slope of the leaf and the 
y-axis is perpendicular to da. That is, n = s . In this frame 


A 

n*k = sin 0 sin q) 


(119) 


Customarily , the measured bidirectional reflection and transmission functions 
are given in the upper leaf frame with h = e* and s down the central line of sym- 

Z X 

metry . When the central line of symmetry is the line of maximum slope , the following 
coordinate transformations are obtained; the quantities that have asterisks are the 
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usual spherical polar anples in the upper leaf frame . 



sin 0 sin cp = cos 0* 

sin 0 cos cp = sin 0* cos cp* 

cos 0 = sin 0* sin cp* 


0* = cos ^ 

sin 0 sin <pj 

u, -1 

r I 

cp* = cot 

tan 0 cos (p 

rv -1 

r . . 

0 = cos 

Sin 0* sin cp 

-1 

r 

cp = cot 

tan 0* cos cp 



(120) 


( 121 ) 


( 122 ) 


In the integration frame , space may be partitioned as shown in table V and 
figure 1 . All L_ radiance incident on the upper surface of da in this partitioning 

is contained in region a , all L_ radiance incident on the lower surface of da is 
contained in region b , all radiance incident on the upper surface is contained 

in region c , and all radiance incident on the lower surface is contained in region 

d . Further , all upwelling intensity emerging from the upper surface is contained 
in region a, all downwelling intensity emerging from the upper surface is contained 
in region c , all upwelling intensity emerging from the lower surface is contained 
in region b , and all downwelling intensity emerging from the lower surface is con- 
tained in region d. (Emerging intensity is that either reflected from the side or 
transmitted from the opposite side , ) These regions give the integration frame limits 
whereby the scattered intensity can be written explicitly gind , hence , the upwelling 
and downwelling scattered power can be obtained, as shown in the calculations that 
follow . The scattered intensities resulting from the downwelling diffuse intensity 
L_ are given from equations (116) and (117) by using the integration limits a and 
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r 


b; thus 


r 

Ipu(k) = -^ypud,h')|n-k'|dfi' 


^ y'ty(k,k')ln.k'|dJ2’ 


b 


L_ r _ 

I^L<k) = — y Y(k,k')|n-k'ldf2’ 
b 


(123) 


(124) 


(125) 


(126) 


The subscripts p and i indicate reflected and transmitted intensity , respectively . 
Similarly , L and U indicate lower and upper surface incidence , respectively , 

Let 1^ denote that portion of Ip^, I^j^, l^^.and directed upward , and 
let 1_ denote that portion directed downward . Then , from the defining equation 
for intensity , the upward-directed power p^ is given by 

2it 

p^ = daj' I^(k) I n*k|dfl (127) 

0 

and the downward-directed power p_ is given by 


2it 

p_ = da J I_(k)|n*k|dn (128) 

0 


45 


Thus , and p_ can be written as 


P+(®n’‘Pn) = 

d. 






A . A 


+ da / |Ipj^(k) |n*k| + I^^(k)|n*k| 


d^l 


(129) 


and 




dO, 




A A ^ . 


k| + I^yCk) |n*k| 


da (130) 


The quantities p^ and p_ are the upwelling and downwelling power from 

a single scattering component in AV . This component is characterized by a slope 

0 . The ensemble of elements in AV has some distribution of slopes . Thus , a 
n 

function S(0^,r|) may be defined so that S(0^,ri)d0^ is the fraction of components 
in AV having slopes between 0^ and 0^ + 6^^ (0 < 0^ < tt/ 2) . The quantity 
S(0n>ri) is a probability density function that describes the distribution of slopes 
and azimuths . Thus , average power <P^> and <p_> are defined by 


<P+> = 



<P_> = 



(131) 


(132) 
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where S(0^,r)) is normalized by 


1 = 



(133) 


Note that and p_ are not irradiances (power per unit area) evaluated on 

a horizontal plane but power per component element directed upward or downward. 
To obtain the power per unit area across a horizontal plane , calculation of the number 
of elements per unit of horizontal area is necessary. A volume element is 


dV = dz dA (134) 


where dA is an element of horizontal area. Then, the number of elements in dV 
is 


N = n dV = n dz dA (135) 


Then, the upwelling power from dV is simply 


dp^ = N<p^> = n dz dA<p_j> (136) 


Similarly, the downwelling power from dV is 


dp_ = N<p_> = n dz dA<p_> (137) 

However, dp^/dA is simply the irradiance on a horizontal plane resulting from 
scattering upward by the collection of elements in dV, and dp_/dA is, similarly, 
the irradiance on a horizontal plane resulting from scattering downward by the com- 
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I 


ponents in dV . Thus 


dp+ 

dE^ = = n<p^>dz (138) 

dp_ 

dE_ = = n<p >dz (139) 


From equations (138) and (139), the change in irradiance on a horizontal plane 
resulting from scattering of downwelling L by the components in the layer dz 
can be written 



= n<p^> 


(140) 


dE 

= n<p > (141) 


Similarly, the scattering of the upwelling uniformly diffuse field can be obtained. 
Let the primes indicate quantities arising from L^. Then 


“ IT ^Pu(k,k')n*k’ df^' 

c 


(142) 


L ^ 

i;b(k) = IT J Tj^(k,k')n-k’ dfi' 


(143) 


" irj pL(k,k’)n*k’ dfl’ 
d 


(144) 


r 

i;u(k) = — yiy(k,k')n*k' dfi' 


(145) 
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The limits c and d were given in table V . 

Dividing the radiation field into upward and downward power flow gives 


p' = da 


/[■: 


py(k)|n.kl + I'^Cklln-klldn 


/t; 


^ , A 


+ da /|rL(k)lh*k| + 1^^ (k) |h*kl | dfl 
b 


|n-klj 


( 146 ) 


/[■p 




p:. = da / |ry(k)ln-k| + l^L(k)ln*:k||dfl 


I n*:klj 


/[■p 


+ da / |rj^(k)|ri*kl + l^y(k)|n-kndf2 
d “ 


( 147 ) 


The quantities <p^> and <p^> are then given from the slope density func- 
tions by 


<Pl> = 



<P> 



( 148 ) 


( 149 ) 
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Thus 


dE’ 

-di- = «<P;> (150) 

and 

dE;^ 

= n<p> (151) 

Thus , change in irradiance , both downwelling and upwelling , due to scattered 
and L_ has been obtained . The next step in calculating expansion coefficient 

ay for i = 1,2 and j = 1,2 is to calculate the total power incident on da from 

L , and L . From this quantity , the total power intercepted by both sides of the 

components in the volume element can be calculated . 

The defining equation for intensity or radiance gives the relationship between 
intensity and power delivered to an area, 

Pj = /ij (6 ,cp)k»da dfl (152) 


A 

where k is a unit vector in the direction of the intensity . In the coordinate frame 
shown in table V 


n»k = sin 0 sin cp (119) 


where n is a unit vector normal to da. When power is intercepted from L_ by da, 

note that L_ is uniformly diffuse and, thus, constant; however, it is not upwelling. 

Therefore, it is nonzero only for regions a an^ b in table V. That intensity contained 
in region a is incident on the upper_side of da, and that intensity contained in region 
b is inciden^on the lower side of da. Thus, the total power delivered to the upper 
surface of da by L_ is given by 
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(153) 


Pu = 


./si 


sin"^ 0 sin q) d0 dcp 


= -d 
2 '' 


+ cos 0) 


and the total power delivered by L_ to the lower surface of da is given by 


p = daL 

1j 


/ 


2 TC 

sin 0 sin cp d0 dcp = ^(1 


cos 0) 


(154) 


Similarly, using regions c and d, the power delivered to the upper surface by the 
upwelling is given by 


P[j = daL^ 


/ 


sin 0 sin cp d0 


dcp = ^(1 


cos 0) 


(155) 


and the total power delivered to the lower surface of da by is given by 


P 


L 


/ 2 7T 

sin 0 sin cp d0 dcp = cos 0) 

d 


(156) 


Adding the power intercepted by both sides , the total power from L_ inter- 
cepted by da is obtained . 


P-P = Pu + Pl 


(157) 


The total power from intercepted by da is 


Pjp = Pu + Pl = da 


(158) 


A peculiarity of a uniformly diffuse field is that the total power intercepted by an 
area is independent of slope or orientation. With equation (158) , the derivation has 
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been completed, and it is possible to write expressions for ^ 21 ’ 

® 22 - 

Returning to the Kubelka-Munk equations and their development , note that 
a^^^E^ is the portion of lost over dz because of either absorption or redirection 

downward of E^. The total upwelling power intercepted by da is given in equa- 
tion (158) . The amount of power emerging upward from da is given in equation 
(148). The difference is that power lost by intercepting da. Thus j 

(159) 

Similarly, ^ 22 ^- change in E_ lost over dz because of either absorption 

or redirection of E_ . Using a similar argument for the derivation of and 

using equations (132) and (157) , &22^- given by 

(160) 


^ 22 ^ 


_ = n L_7 c 


da - <p 


^ 11 ^+ 


= n L^tt da - <p^> 


The expression portion added to E^ over dz by scattering of E 

into E^. This quantity is given directly from equation (140) by 


dE^ 

®12®- ' -3T = 

Similarly, is the change in E_ over dz resulting from scattering of E 

into E_ . This quantity is given directly from equation (151) by 

^21^+ ~ (162) 


i 
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Note that E = h % and E. = L it. The first four a., can be written 
- - + + 1] 


^11 ” 


da - 


<P+> 


(163) 


n<p^> 
‘l2 ^ L 7t 


(164) 


n<p> 

^21 ^ L^Tt 


(165) 


®22 ^ 


da - 


<P_> 

L 7C 


(166) 


Equality of some of the a^ is obtained 


The preceding a„ are not necessarily equal. 

only under speciM circumstances to be discussed in a later special-example calcula- 

a^g undeveloped. These 


tion. The preceding derivation leaves only a^^g. 


and 


coefficients relate changes in the specular or directed component incident on the 
upper surface of the canopy . Dealing with the specular component is somewhat more 
complicated because the incident intensity is not uniform over a hemisphere; however, 
it is simplified somewhat because integrals over beams (differential solid angles) 
are trivial. 


Development of Expansion Coefficient 

A primary assumption for all models has been that the scattering components 
(leaves , branches , etc.) are representable by a collection of planar scattering ele- 
ments da, each of which has an associated slope 0^ and azimuth cp^ and a unit 

normal to the surface n . Then , the slope angle 0^ will be the customary zenith 

angle 0 in a z-axis vertical coordinate frame . Furthermore , density functions may 
be defined as D(cp^,v) (eq. (42)) and S(0^,q) (eq. (43)). The expression 

D(cp ,v) represents the fraction of leaves having an azimuth between cp^ and 

cp^ + expression S(0^,ri) represents the fraction of leaves having a 

zenith angle between 0^ and Using the outward- drawn unit normal, 

a vector element was defined in equation (5) . 
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Define Iq(0q,cPq) as the specular intensity incident on the upper surface 

of the canopy; this intensity will be contained in the solid angle intercepted 

by the Sun. The expression ^^(9^,9^) is a unit vector in the direction of the solar 

or specular intensity. Then, the irradiance on a horizontal plane at the upper surface 
of the canopy would be given by 


E^(0) = (167) 

where e cos 0 q and where the irradiance of the specular beam measured 
perpendicular to the beam is given by 


^b = ^0 ^"^0 


(168) 


The cross section presented by some da to the beam will be given by 


da'kg = da [cos 


(169) 


The expression for cos y^ is given by considering the dot product of and 

n expressed in a z-axis vertical coordinate frame where the x direction J^direction 
of 9 = 0) is left arbitrary. (See eq. (41) .) Defining the angle between k^ and 
n to be y^ gives 


I cos y„| = I sin 0 „ sin 0 cos l(p^ - 9 1 + cos 0 „ cos 0 I (170) 

' 'O' ' 0 n (^0 ^nJ 0 n' 

The expression for cos y^ is then a function of 0^ and 9 ^ explicitly. Using the 

density functions D (9 ,v) and S(0 ,r|) , an average cross section may be defined 
by n n 


IL 

2n 2 


<da> 


9 0 
^n n 




COS y. 


|d 9 d 0 
' n 


(171) 
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or 


0 "" da<|cos Yq|> 

n 


(172) 


where 


— 

2Tt 2 

<|cos YqI> = f f D^9n.v^S^0j^,Ti^|cos (173) 

where it is assumed that no leaf has a slope >7t/2; that is, elements do not turn over. 

Performing the integration over 0^ and cp^ to obtain <cos Yq> is completely 

straightforward if cos Yq has no negative values . If cos Yq has a negative value , 

as would usually be expected, and thus passes through zero, the integral must be 
broken up into a sum of integrals over regions for which the sign of cos Yq does 
not change. 

The average value of a function i (t) over an interval T is given by 

T 

^ f ()(t)dt (174) 

0 


As noted , < | f) (t) | > = < fj (t)> if and only if the sign of j$ (t) does not change . 

If, however, there are regions of the interval T for which ^ (t) possesses a nega- 
tive sign , the absolute value of the function over the whole interval T is found by 
taking the negative value of the integral over the regions in which the function is 
negative and the positive value of the integral over the regions in which the function 
is positive . 



J l$(t)dt - J Ut)dt 


T+ T_ 


(175) 
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where denotes the positive regions of ^(t) and T_ represents the negative 

regions of |j(t) within the interval T . The boundaries of these regions are marked 
by the function passing through zero; thus, the boundary is found by solving 
/j(t) = 0 to obtain Tq , 

For the specific instance of obtaining <da>, the basic problem is averaging 
cosYq, Setting cos Yq = 0 (eq- (170)) gives 


sin 0 sin 0 cos/^cp - cp '\ + cos 0 cos 0 = 0 
n o l^^n ^oj n o 


(176) 


or 


cos 6 = -cot 0 cot 0 (177) 

n o 

where 6 is defined as the value of 9^^ ~ for which cos Yq = 0. Physically, 

this value corresponds to the value of 6 for which the leaf is parallel to the solar 

beam. Because 0 and 0 are O<0 <%/2 and O<0 <-kI2, respectively, 

n o -n~ -o- ^ 

both cot 0^ and cot 0^ are positive . Thus , because the cosine is never greater 
than 1, the condition for 8 to exist can be expressed as 

cot 0 cot 0 <1 (178) 

n o - 

from which the following may be written 

|<0nH-0o<^ (179) 


as can be seen from the following trigonometric manipulations . 


cot 0 cot 0 
n o 


cos^0 - 0 ) + cos(0 + 0 ^ 

V n o) \ n o) . . 

cos/0 - 0 ) - cos/0 + 0 ) - 

V n o) V n o; 


(180a) 
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or 


cos 


(®n - ®o) * ®o) ^ “=(®n ' ®o) “ * ®o) 


or 


cos 


(®n * ®o) 5 


(180c) 


Thus, 


? < 0 +0 < 7t 

2 - n o 


(179) 


The average of cos Yn over cp can then be written, using equations (178) and 
(179), as u n 


<|C0SY„I> =25- 

^n 


8+cp 27t-6-cp 2n 

f °-/ 

o 6+cp 27t-6-cp 


ol 


cos Yq ^^'Pn)^(^n) 


(181) 


Then , averaging over 0^ is done by rewriting the existence condition for 6 (as 

defined in eq. (179)) as 8^0 if and only if 7t/2 - 0^ < 0^ < tt/ 2 and as 6 = 0 

if and only if 0 < 0^ < n/2 - 0^. Thus, two ranges over 0^ are established; 

that for which 6 exists and that for which it does not. Then, the complete expres- 
sion for the average of cos Yq over 0^ and 9^ is given by 


1 7t_„ 7t 

1 2 '^o 2Tt 2 

If r r 


\ 


27t-6-cp 2k j 

r ° r 1 

J J " / 

f - 

/ "i 

(0 0 

1 2 ®o 

J 

0 

8+cp 27t-8-cp 1 

^0 9 JI 




cos 




(182) 


57 



For the interesting phenomenon of azimuthal symmetry where D (9 ) 
equation (182) can be simplified by using the general relationship ^ 


2it 


/ 


cos vl/ d\|/ = 0 


and the normalization equation for the azimuthal density function (eq. (44)) 
forming the integrations gives 


<l°°®Yol>e cp =4 
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However, equation (177) and the identity tan = sin/cos give 


cos 0 = 

n 


sin 0 cos m + 9 1 sin 0 
n \ ^o) q 

cos 0 


1 , 


(183) 


Per- 
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Thus , the second and third terms of equation (184) can be combined . 


q> 


n’^n 



cos 0 d0 
n n 



(186) 


The quantity a^^ can now be developed. The amount of Iq(0q,cPq) inter- 
cepted over a path d2 along will be proportional to the fraction of the 

total cross section of the beam obscured by the elements in a volume element along 
dfi . A volume element dV along dC is given by 


dV = da de 


(187) 


where da is the cross section of the volume element . The number of elements in 
dV is given by 


N = n dV 


(188) 


where n is the number density of elements . The total cross section of all elements 
in dV is then 


N<da> = n dV<da> 

However, the fraction of the beam obscured is given by 

F = = n dC<da> 

da 

Then 


dE 

dC 



= n<da>EQ 


(189) 


(190) 


(191) 
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However 



(192) 


Therefore 


n<da> 



( 193 ) 


Now , only the quantities a^^^ and must be developed . This development 

is complicated for several reasons . The specular incident intensity is incident on 
the top of some leaves and on the bottom of others. Because 1^ depends on cp , 

azimuthal dependence becomes important . For this reason , all elements , including 
the bidirectional reflection and transmission functions , must be referenced to a sin- 
gle coordinate frame. 


Development of Expansion Coefficients 



and 


The expansion coefficients a^^^ and ^ 23 ’ relate the specular intensity 

to that scattered upward and downward, are obtained by noting £^^^ 3^0 as the change 
in upwelling diffuse intensity over dz due to scattering of specular intensity where 


^0 *0 ^^ 0^0 


( 194 ) 


Similarly, a23®0 downwelling diffuse intensity over dz due to 

scattering of specular intensity . 

From the development of the phase function , it was shown that the scattered 
intensity from a volume element upon which intensity Iq was incident was given 
by equation ( 22 ) . 


Is(k) 



An, 


( 195 ) 
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From equation (15) , the power (p^ for upwelling and p_ for downwelling) scat- 
tered from the volume element is given by 

p^ = AA f Ig da = j- AVIq AHq J P^k.kojda (196) 

where ^h refers to upper or lower hemisphere, respectively . Thus, p^ is up- 
welling and p_ is downwelling power . However , 


P± 

AV 



(197) 


and 


dE^ 

dz 



(198) 


Therefore , a^^^ and are given by 


1 

^13 



^ 1 

^23 47tpQ 



(199) 


( 200 ) 
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Bidirectional Reflection Function for the Canopy as a Whole 

The objective of the Kubelka-Munk approach is to obtain an estimate of the 
radiation field at any level within the canopy . This estimate is then used in the 
formal solution of the radiation transfer equation to obtain the upwelling field at the 
upper surface of the canopy . The formal solution to the radiation transfer equation 
as shown in the appendix equations (269) to (278) may be written 


P dl(t,k) 
dt 


I(t,k) 


F(t,k) 


( 201 ) 


where i is the optical depth of the canopy, with F(x,k) given by 


1 2n 

F(t,k) = ^ JJ P(k,k')l^(t,k')di^’ 

-1 0 

The formal solution is given by the method of integrating factors as 


( 202 ) 


‘0 I 

where is the upper surface of the canopy and x' is specific optical depth , with 
dt = -n <da> dz and k_^ being the upwelling or the downwelling k. The function 
l^(t,k') is the total field at c, which is obtained from the Kubelka-Munk solution. 
Thus 





T-t' 




dt' 


(203) 


(t . k) 


L^(t) 


+ L (t) + 



(204) 
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However 


L^(t) 


E^(t) 

% 


L (X) 


E_(x) 

7t 


and 



Thus, F(x,k) can be explicitly written 



where 


Xq is the upper surface of the canopy . 


(205a) 

(205b) 


(206) 


(207) 
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From a knowledge of I (Tq , k) , it is possible to specify a bidirectional reflec- 
tion function for the canopy as a whole . The defining equation for the canopy bi- 
directional reflection function can be written 


Sy(T;k,k') 


47I|I 


l(0,k) 


(208) 


where l(0,k) is the upwelling intensity at the top of the canopy, p is the cosine 
of the zenith angle of Ic , and t is the optical thickness of the canopy as defined in 
equation (275) in the appendix. The expression I (k„) represents the illumination 

at the top of the canopy contained in The canopy bidirectional scattering 

functions are discussed in greater detail in appendix equation (288) and those that 
follow . 

Once the canopy bidirectional reflection function is known, the upwelling ra- 
diation for any condition of illumination may be calculated . Because the canopy bi- 
directional reflection function is so complicated , a thorough understanding of its 
characteristics can only be attained by programing the model and doing variational 
studies . However , valuable insight can be obtained by considering some less com- 
plex , special examples . 


MODEL II 


The purpose of the canopy model is to aid in the interpretation of spectral 
data. Ideally, the probably unattainable objective for a canopy model would be to 
provide unique , theoretical signature inverses from which signatures or collections 
of signatures could be clearly interpreted in terms of type , condition , maturity , 
yield, and so forth. Realistically, however, the model should improve understand- 
ing about how various characteristics affect signatures or how data might best be 
taken to minimize ambiguities . In short , a realistic goal would be that the model 
provide understanding and , thus , improve data-taking and analysis techniques . 

One problem in using the previously derived model approach is the complexity . 
It is difficult to gain insight simply from the study of the structure of such compli- 
cated systems of equations . Insight can be gained through programing the models 
and performing parametric variation analysis . 

Another practical problem is the paucity of bidirectional scattering data for 
canopy components . Although some information is available in the NASA Earth Re- 
sources Spectral Information System , it is not sufficient . The best available work 
known to the author was done by Breece and Holmes (ref. 3) , who have studied the 
bidirectional scattering characteristics of green soybean and corn leaves for top 
incidence . 
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The obvious approach is to consider some simpler but , hopefully , sufficiently 
realistic models . At best , the simpler models based on reasonable assumptions (at 
least for some wavelengths) may produce solid predictions of canopy reflection . 

At worst , use of such models will simplify the equations so that insights for guiding 
the parametric equations may be gained by merely studying the structure of the 
simpler equations. 

The most apparent simplifying assumption is that the bidirectional scattering 
functions are constant; that is , to assume that the components are uniformly diffuse 
or Lambertian scatterers . Breece and Holmes suggest that this is an excellent 
assumption for the transmission function at wavelengths from 0.37 to 1 micrometer 
and almost equally as good for the upper surface reflection function at wavelengths 
higher than 0.75 micrometer. An array of planar scatterers, each of which is indi- 
vidually Lambertian , does not necessarily give rise to Lambertian scattering for 
the entire array. 

Additional simplifying assumptions for the modified Chandrasekhar approach 
include assuming the azimuthal density function D(cp^,v) to be constant and the 

slope density function S(0^,ri) to have some analytic form, such as a Gaussian 

distribution around some mean slope . A further simplification is obtained by con- 
sidering only the vertical look direction, which is appropriate to spacecraft-mounted 
sensor systems . 

The first special situation is that of constant leaf bidirectional reflection func- 
tions . The reflected and transmitted intensities for upper surface incidence on com- 
ponent da are found from integrating equations (123), (124), (142), and (145); 
for lower surface incidence, from integrating equations (125) , (126) , (143) , and 
(144). In the following equations , G is defined as (7t/2) (1 - cos 0) and F is 
defined as sin 0 sin cp . 



(209) 


y r dn 

b 


= G 


( 210 ) 



dfl 


= G 


( 211 ) 


/ 


r dfi 


7t - G 


( 212 ) 
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where 


dfl = sin 0 d0 dcp 

Thus, for L_ , integrating equations (123) to (126) gives 



L_P 


U 


% 


1% - G] 


(213) 


(214) 


‘,L = CG] 


Vl 71 


L-'tt 

'tu = -ir - G] 


Integrating equations (142) to (145) for gives 


r = 

pU Tt ^ 


L^r 


tL 


*(.-«) 


^pL 
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(215) 


(216) 


(217) 


(218) 


(219) 


( 220 ) 
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These expressions can then be used to integrate equations (129) and (130) 
to give upwelling power and downwelling power due to downwelling incident 
intensity . 


P+ = ^ - G) + tiP] Or - G) + j^PjG + - G)Jg| 

P_ = ^ - G) + Pjp\ (^ - G) + [tjG + PyOr - G)]g| 


( 222 ) 


(223) 


Integrating equations (146) and (147) gives upwelling and downwelling power due 
to upwelling incident intensity 


Pi = 


|[tL(Tr - G) + p^Gj (7t - G) + [tyG + Pj^(7r - G)]g j 


Pi = + x^g] (X - G) + [pyG + ij^iTZ - G)]g| 


(224) 


(225) 


Then set 


TtD' 



7tR' 



(226) 


(227) 
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Thus, using equations (131) , (132) , (148) , and (149) , the average upwelling and 
downwelling power from scattering of downwelling and upwelling intensity is given 
by 


<P+> “ baL 


-|[Pu Pl " - PuJttR' + PyTcj 


(228) 


<p_> = daL_ 


2t - ^Pu + PlJ + Pl Pu ‘ 


<p]^> = daL^ 


j “ (Pu Pl^ I^Pl Pu ~ + T7t| 


<p^> = daL 


j[pu ^ Pl ■ 2']’'“’ - Pl)] Pl”| 


(230) 


(231) 


Then, using equations (163) to (166) , which define the expansion coefficients, 
^11’ ^12’ ^21’ ^22 be written 


a , = n da)l - t 

'' { 


2t 


■12 


21 


22 


(pu ^ Pl)]*’’ * [(pu * Pl) ■ 

n daj^l^p^ + - 2,jD’ + 2^, - p^] R' + p^j 

n dajpj^ + + pj^ - 2 i1d' + 2j^t - PjlR'j 


n dajl - T + j2i - ^Py + Pl^J jj^PL Pu^ " 2^R'^ (235) 


(232) 
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Further, consider the situation in which Py = pj^; this gives 


11 


= n dA|l - T + (t - p)b| 


( 236 ) 


a^2 = dA|p + (p - t)b| 


(237) 


= n dA|p + (p - t)b|. 


(238) 


22 


= n dA^l - 


T + (t 


P)b| 


(239) 


where 


D 


t _ 


7T 


R' = B = 


l/°(l-)s(vn) 


d0. 


n 


(240) 


Notice that where p = t , the B expression , which carries all the slope characteris- 
tics , is removed from the equation . A preliminary conclusion is that in those spec- 
tral regions in which transmission is nearly equal to reflection , the effects of leaf 
slope are minimized; thus , variations in slope caused by wind should cause decreased 
statistical variance in S(k,k') and , hence , in signature . 

The expansion coefficients a^^ and 32^ for constant bidirectional scattering 

functions can be found in a straightforward manner by noting that the reflection and 
transmission scattered intensity is proportional to the irradiance . Thus , from the 
defining equations for the bidirectional scattering functions (eqs. (37b) and (116) 
to (122) and table V) , the intensity due to reflection from the upper surface ^pu» 

the intensity due to transmission from upper to lower surface 1^^ , the intensity due 

to reflection from the lower surface , and the intensity due to transmission from 
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the lower to the upper surface can be written 


pu 
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pL 
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where 


k(,.n i 


COS Yq = sin sin 
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0 cos 
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(<Po - %) 


cos 0Q cos 


0 


n 


(245) 


If cos Yq is positive, 
of da and = !,„ = “• 


then incident intensity is on the upper surface of da 
is negative , then incidence is on the lower surface 


In the section entitled "Canopy Geometric Density Functions," it was shown 
that the conditions for upper and lower incidence can be written as follows . When 
cos Yq > 0 


0 < 0^ < I 

n ^ 


‘PQ 


0 < cp < 2n 


(246a) 


or 



< 0 


n 





6. < cp < cp. + 8. 
1 ^n ^1 1 


(246b) 
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When cos < 0 


7t 

2 


e„ < e 

0 n 



CD. + 6. < q) < cp_ + 2 ti - 6. 
^1 1 ^0 1 


( 247 ) 


where , expressed in a z-vertical coordinate system 


6. = cos cot 0^ cot 0g (248) 


The upwelling power reflected from the upper side of da is given from a considera- 
tion of figure 1 by 


PpU “ Ipu da ^ cos Y dfi (249) 

a 

where the upwelling power is independent of coordinate representation; therefore, 
a convenient representation is chosen: n = x and y is along the slope of da. 

Thus, 


cos Y = sin 0 sin <p (250) 

Similarly , the downwelling reflected power from the upper surface is given by 

Ppu = ipu y* y 

c 

The downwelling transmitted power is given by 

p^U “ ^tU J' ^ (252) 

d 


and the upwelling transmitted power is given by 

Ptu ^ *tu f y 

b 


( 253 ) 
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Combining equations (249) , (253) , (241) , and (244) gives the total upwelling 
power for a scattering element da illuminated on the upper surface. 


Pu 



‘o ''"o 

cos Yn da 

7C ' 0 


Pu /■" "u / 1''°® Y dn 


(254) 


And combining equations (241) , (244) , (251) , and (252) gives the total downwelling 
power for a scattering element da illuminated on the upper surface. 
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/ \ ^0 ^^0 1 
(0 ,9 1 - " "da cos Yn 

\ n ^nj 7c '0 



'' ' 1 

c 

d 


cos Y dn 


(255) 


Note that these expressions are for a particular da characterized by particular 

values of 0 and cp . 

n ^n 

The integrals appearing in the preceding expressions have been evaluated 
in equations (209) to (212) and, thus 


+ 




cos Yq da 
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and 
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Iq ^^0 ''^0 
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(257) 


where, by definition, G = (tt/ 2) (1 - cos 0) . When elements are illuminated by 1^ 

in AHq incident on the lower side of the element da , they are handled similarly . 

The reflected intensity is given in equation (243) , and the transmitted intensity is 
given in equation (242) . Then, the upwelling reflected power is 
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pL 



da 


J" kj.*n dfl 
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and the downwelling reflected power is 


p ^ = I T da / k *n 

^pL ph J r 




(259) 


The upwelling transmitted power is 


PtL ' ‘tL 


/ A A« 

k^.n 


dfi 


(260) 


and the downwelling transmitted power is 


p^=I_da fk-n 
^tL tL J ^ 


dJ2 


(261) 


Using the already evaluated integrals in equations (209) to (212) , the total upwelling 
power for lower incidence from da is 




lo "^^0 

_ 


da 




(262) 


and the total downwelling power for lower incidence is 



In ^^n cos y. 


TZ 


da 


pL” * pL - Pl) °] 


(263) 
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The total upward-directed power is obtained by integrating all da over the 
distributions expressed in the density functions S(0^,r|) and D(cp^,v) , where 

the limits are obtained from the conditions on the sign of cos y- given in equations 
(246) and (247) . " 

The total upward-directed power by all elements in AV is given by 


p = n AV 
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and the total downwelling power by all elements in AV is given by 
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As in equations (196) to (200) , and a^g found from the total upward 

and downward power by 


a = -B- 
®13 AVE, 


(266) 


1 = 

23 AVE, 


(267) 


Because the quantity a^g is independent of the bidirectional scattering functions, 

it is unchanged from the expression given in equation (193) . Obtaining the same 
upper and lower bidirectional reflection functions simply requires removing the sub- 
scripts U and L from the previous expressions for expansion coefficients a^ . 

Interestingly, in the wavelength regions where p = t, the geometrical factor 
G disappears . Thus , the effect of leaf slope in those regions may be expected to 
be diminished inasmuch as leaf slope appears only in Sgg. Therefore, in these 

spectral regions , leaf slope variations probably will be dominantly seen in the single 
scattering component of the specular intensity as expressed in the agg term. In 

addition, slope will be expected to affect optical depth as defined in equation (275) 
(appendix) , which will then change the effect of lower layers and the soil. However, 
if the canopy optical depth is a factor of 3 or more (that is , the specular beam has 

_3 

been reduced by a factor of e ) , then slope will have little effect in exposing the 
soil to either illumination or view . Thus , soil will have little effect on reflected 
intensity . 


MULTILAYER CANOPIES 


Heretofore , the canopy model calculations were explicitly written for a single- 
layer, single-constituent model. For some canopies, these calculations will be suf- 
ficient . For others , extension of the model to multilayer , multiconstituent canopies 
will be necessary. All the fundamental calculations and principles for the multilayer, 
multiconstituent canopy were developed in the single-layer, single- constituent cal- 
culations. An explicit expression of the multilayer, multiconstituent model would 
require an extremely complicated notation. However, a description of procedures 
for performing the extension can be easily presented and should be sufficient for 
any investigator attempting to program the model . 

Generalization to multilayers affects the model at essentially two points: the 
boundary conditions on the Kubelka-Munk equations for a given layer and the calcu- 
lation (using the formal solution to the radiation transfer equation) to find the up- 
welling intensity at the surface of the canopy . The boundary conditions for a single- 
layer canopy are used to solve for the terms Q and R as given in equations (95) 
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to (98) . In the multilayer calculation, each layer is subject to a separate set of 
Kubelka-Munk equations , constructed for each layer in exactly the same way as for 
the single-layer canopy . Thus , each layer will have Q and R terms . The set 
of equations resulting from the equations for Q and R for each layer will be a set 
of simultaneous algebraic equations that must be solved using the boundary conditions 
appropriate for the multilayer calculation . The multilayer calculation gives a set 
of equations in Q. and R^ where i is the index for the layer . If there are four 

layers , then i takes the values 1 to 4 . 

The boundary conditions at the top of the canopy and at the canopy-soil inter- 
face are the same for one layer or for many layers . These boundary conditions are 
given in equations (101) and (102) . An additional boundary condition between layers 
is that the downwelling uniformly diffuse radiation at the bottom of one layer is equal 
to the downwelling uniformly diffuse radiation at the top of the layer immediately 
below it . Similarly , the upwelling uniformly diffuse radiation at the top of any one 
layer is equal to the upwelling uniformly diffuse radiation at the bottom of the layer 
immediately above it. The attenuated or reduced specular beam intensity at the bot- 
tom of one layer is equal to the attenuated specular intensity at the top of the layer 
immediately below it as well . In short , the total radiation field at any interface be- 
tween layers has only one value , and that value is shared by both layers . 

The other point at which multilayers affect the model calculation is in the use 
of the formal solution to the radiation transfer equation (eq. (2)) to find the upwelling 
intensity at the top of the canopy . The problem may be most simply solved in a 
sequential, layer by layer manner, beginning at the bottom and working upward. 
First , the radiation transfer equation formal solution is used to find the intensity 
at the top of the bottom layer by using the Kubelka-Munk-derived fields for the bot- 
tom layer in precisely the same way described for the single-layer model. By impli- 
cation, a phase function for the first layer must be calculated by using precisely the 
same procedure described for the single-layer model. Once the intensity at the top 
of the bottom layer is found, this value becomes the same as IqC^) for the second 

layer, and the process is repeated for each successive layer until the value of up- 
welling intensity at the top of the top layer is found . This value of upwelling inten- 
sity in the k direction is then used in precisely the same way to determine the total 
canopy reflection scattering function S(k,k') as was done in the single-layer 
development . 

Generalization to multiconstituent canopies affects the model development in 
two places . The first place is in the a^ , the expansion coefficients for the Kubelka- 

Munk equations . The second place is in the development of the phase function . 
Conceptually , the generalization of the a.^ to multiconstituents is simple . Because 

the a., are developed from the assumption of single scattering within some volume 

element , a separate a. . can be developed for each constituent , These a., are then 

1 ] ^ 

added. If the constituent index is taken as m, then for an N-constituent layer, the 
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total a. . 
1 ] 


for the layer would be given by 



( 268 ) 


The total 


a.. 

1 ] 


would then be used to solve the Kubelka-Munk equations in the same 


way as for a single-constituent layer . Similarly , the phase function for a multicon- 
stituent layer would be the simple sum of phase functions for each constituent . 


Clearly , conceptual generalization of the model to multilayers and multicon- 
stituents is not particularly difficult; however , the notation and computation can 
become complex. In summary, computation for multilayers and multiconstituents 
within layers requires the calculation and summing of a set of a.^ and the calcula- 
tion and summing of a phase function for each constituent within each layer . Of 
course, the physical system could enable some fortuitous simplifications. For exam- 
ple , if the canopy consisted of leaves in which only the volume density varied with 
height and if reflectivity , transmissivity , and slope and azimuth distributions re- 
mained the same , the calculation of the a_ and the phase function for multilayers 

would not be difficult . Similar simplifications might be effected if only slope or 
azimuth distributions changed with height or if only reflectivity and transmissivity 
changed with canopy height . 


UNRESOLVED QUESTIONS AND DIRECTIONS FOR FUTURE WORK 


Although much has been done , much remains to be done in the fundamental 
theoretical analysis of canopy reflection . A differential equation approach requires 
the assumption that a differential volume element can be meaningfully defined. This 
assumption implies the existence of a volume element that is small compared to the 
depth of the canopy and that has average values similar to those of adjacent volume 
elements . The validity of this assumption seems to be related to leaf density in a 
way that is not quantitatively clear. Thus, one important question is over what di- 
mension can homogeneity be assumed? A related question is what canopy character- 
istics are necessary for the Kubelka-Munk-type approximation to be valid? These 
characteristics are probably a function of both density and geometric distribution . 
Such questions could be usefully investigated in a laboratory using a carefully con- 
structed artificial canopy having translucent plastic arrays . 

Another unanswered question concerns the horizontal inhomogeneity or row 
effects on canopy reflectance . Purely theoretical consideration of this problem is 
even more difficult than that for the nonrow example considered in this paper. 
Hopefully, the study of special, less complex situations will enable simplification of 
the problem to the point at which the horizontally anisotropic situation can be made 
tractable . 
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Another important and possibly critical question that is not addressed in this 
paper or in other current modeling efforts is the theoretical characterization of grain 
heads such as those found in wheat . The straightforward characterization of a head 
of wheat as a collection of planar elements may not be possible . This question is 
important because 50 percent of the top layer volume of mature wheat can be occupied 
by heads. The solution to this problem may be simpler than it appears; however, 
the problem remains unsolved. 

In the future , a workable canopy model less complex than those presented 
in this paper probably can be developed . This development should be possible once 
the dominant effects are clearly established . For example , it may be found that the 
radiation reflected only once from the leaves is the dominant effect and that the re- 
mainder of the upwelling radiation of the canopy can be ignored . Further , because 
single-surface radiation often is strongly polarized , it is probable that front surface 
scattering can be isolated and also that polarization may prove to be a powerful remote- 
sensing diagnostic tool . Even though the models could be generalized to include 
polarization by writing them in polarization matrix form , this formulation has not 
been done . Clearly , analysis of the models presented in this paper raises as many 
questions as it answers . 


CONCLUDING REMARKS 


An ultimate usefulness of canopy modeling in remote sensing is the identifica- 
tion and quantification of the causes of variation in the statistical properties of canopy 
signatures. These statistical properties include the mean vector, the covariance 
matrix , and the correlation matrix or some other measure of separability . Excluding 
atmospheric effects , the models presented can be used to identify these causes of 
signature variation. Little can be said quantitatively, however, until models have 
been programed and experimentally verified, factorial tests have been performed, 
and limits on applicability have been determined . Nevertheless , some insight can 
be gained from the structure of the model equations . 

The causes of canopy signature variation are identified as variations in com- 
ponent (leaf, stem, tassel, etc.) slope and azimuth distributions as expressed in 
the slope and azimuth density functions , variations in component density , variations 
in component bidirectional reflection and transmission , variations in soil reflection 
properties, and variations in illumination conditions of both solar and sky components. 
The complexity of the physical situation is evident; thus, the complexity of the canopy 
model follows. However, in spite of the complexity, analysis of the models shows 
that the canopy modeling problem is solvable. Because of the complexity of the 
physical situation , it is probable that a purely experimental program to determine 
the cause and effect relationships of canopy signatures is impractical and may be 
impossible. Thus, the use of theoretical models to quantify signature variation cause 
and effect appears essential in any program. Further, because of the complexity 
of the physical situation expressed in the models, it is likely that experimental pro- 
gram data and remote-sensing measurements will be taken under circumstances that 
result in simplifications . One possible simplification is to restrict measurements 
to that portion of the spectrum in which skylight is minimal and to atmospheric con- 
ditions in which cloud contributions are low . Another possible simplification would 
be to restrict look directions to the vertical , though such a restriction may impose 
undesirable limits on diagnostic capability. 
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A final conclusion from the structure of Model 11 is that there are probably 
regions of the spectrum in which variations in slope and azimuth distributions will 
have minimal effect on signature variation . In these regions , the components are 
approximately Lambertian and the transmission and reflection of the leaves is approx- 
imately equal. 

Although this investigation is incomplete , the insights presented will hopefully 
contribute to raising the quality of future remote-sensing research . 


Lyndon B . Johnson Space Center 

National Aeronautics and Space Administration 
Houston, Texas, March 2, 1975 
951-16-00-00-72 
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APPENDIX 


MODEL III - THE GENERALIZED CHANDRASEKHAR APPROACH TO SOLVING 
THE RADIATION TRANSFER EQUATION FOR A PLANAR CANOPY 

INTRODUCTION 


Canopy modeling must relate incident intensity to upwelling intensity at the 
upper surface of the canopy . This relationship is obtainable from the solution to 
the radiation transfer equation for the canopy . Solving the radiation transfer equa- 
tion has two major aspects. First, a phase function must be developed; second, 
the total intensity field at every point in the canopy must be found . A technique 
for solving the radiation transfer equation that is suitable for canopies and that can 
be approximated by plane-parallel homogeneous layers is discussed in this appen- 
dix . The technique , Model III , is a modification of a technique first proposed by 
Chandrasekhar in his book "Radiative Transfer" (ref. 1) . Therefore , the notation 
used herein will correspond to that used by Chandrasekhar as much as possible and 
only the details necessary to clarify the modification and the reasons for it will be 
included. Complete understanding may require familiarity with Chandrasekhar’s 
"Radiative Transfer" (ref. 1) , in particular with sections 1 to 9, 13, and 49 to 51. 

In the Chandrasekhar technique , the problem is essentially stated in terms 
of new variables , then the radiation transfer equation is used to solve the reparam- 
eterized problem. These new variables are the bidirectional reflection and trans- 
mission scattering functions for the canopy layer as a whole . 

The primary reason for modifying the Chandrasekhar technique is the dis- 
similarity between the canopy and the atmospheric phase functions . Unlike atmos- 
pheric phase functions, the canopy phase function (1) has no symmetry around the 
incident direction , (2) has an orientation with respect to the vertical , and (3) does 
not have mirror symmetry with respect to the horizontal (that is , scattered intensity 
due to radiation incident from a direction bears no necessary symmetric relationship 
to radiation scattered from a beam incident from the opposite direction) . 


RESTATEMENT OF THE RADIATION TRANSFER EQUATION 
FOR PARALLEL LAYERS 


Before restatement of the canopy problem in terms of bidirectional scattering 
functions , it is useful to reexamine the radiation transfer equation and write it in 
a form particularly suited for plane-parallel layers . This form has great similarity 
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to the standard form of the radiation transfer equation for plane-parallel atmospheres . 
The general radiation transfer equation is 


+ n(C)o(C; p,cp)I(e; p,9) = nJ(8;p,cp) 


( 269 ) 


Note that the dependence on p and cp is given explicitly rather than implicitly .by 
using the unit vector The formal solution is obtained by the method of integrating 
factors and shown to be 


I(C; p,9) 




J(E; p,cp) 
a(E; p,9) 


g-T(£,E') 


dE' 


( 270 ) 


where 


T 




n(E)o(E; p, 9 )dE 


(271) 


and J(E;p, 9 ) is given by 


J(E; p,9) 


4ir 


1 2it 

// I^(C; p' ,9')F (E; p,9; p' ,9')df2' 
-1 0 


(272) 


The assumption of homogeneous parallel layers allows change to be expressed 
in terms of the vertical direction z . Thus , 


p dE = dz (273) 

where -1 < p < 1 and where p = cos 0. The variable F is defined by 


F = 


J 

d 


(274) 
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and a differential optical depth is defined by 


dt = -no dz 


(275) 


Then , the radiation transfer equation may be written 


\i dl^(t;p,cp) 
dt 


I^(t;|a,<p) - F(t;n,cp) 


(276) 


where the subscript t is used to emphasize that the expression applies for the total 
field . The variable F is then given by 


1 2 % 


F(t;p,(p) = ^ J j' P(p,cp;p',9’)I^(t;p',9’)d|a' dcp’ (277) 


-1 0 


This form of the radiation transfer equation is formally solvable using the integrating 
-(t-tQ)/ix 

factor e , which gfives the solution 



where the absolute values ensure that regardless of the direction in which the depth 
T is measured, the solution is correct. 

Again , the physical meaning of the solution is easily seen . The upwelling 
intensity at t, which is I(t; +p,cp) , is a result of upwelling intensity at t„, which 

is I(tQ; p.,9) , decreased by e plus that scattered into the beam along 

the path decreased by e ^ ^ where t > i’ and where x - t' is the optical 

distance between the volume element doing the scattering and x integrated (summed) 
over eill volume elements along the path between and x. Although this solution 

will be found useful in the modified Kubelka-Munk approach to canopy modeling , 
it is not in itself a solution for the canopy as a whole. 
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The S and T Functions for a Planar Medium 

Radiation incident on the surface of a layer of some optical thickness will be 
scattered in part , absorbed in part , and unaffected in part by the layer of material . 
The absorbed portion is that removed by any process from the wavelength interval 
of interest and is thus without interest for a reflective model , The scattered portion 
may be directed back into the hemisphere containing the incident radiation , or it 
may emerge from the surface opposite that of incidence . These two scattered portions 
will be designated reflected diffuse radiation and transmitted diffuse radiation, re- 
spectively . (Reflected diffuse radiation is equivalent to what Chandrasekhar called 
scattered diffuse radiation.) In this context, diffuse does not necessarily mean uni- 
formly diffuse, perfectly diffuse, or Lambertian. The radiation that does not encounter 
scattering centers in the layer will emerge from the other side of the layer unchanged 
in direction but reduced in magnitude . This portion of radiation is conveniently 
treated separately and is labeled reduced incident flux. 

Let the total scattered or diffuse portion of the field at t be designated by 

I (t;p, 9 ). Then 

s 

I (t;p,cp) = I (x;p,<p) (279) 

o U 


where 0 < p < 1 , and 

= i^j(t;^.<p) (280) 


where -1 < p < 0 and where is the upwelling diffuse radiation and is the 

downwelling diffuse radiation. Let I.(t; p,<p) represent intensity incident on a layer 
at T. By using the identity 

b -a 

^(x)dx = y* jj(-x)dx (281) 

a -b 
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the defining equations can be written with all functions specified over the interval 
0 < n < 1 to conveniently establish that +p in the argument of a function represents 
radiation in the upward direction and that -p represents radiation in the downward 
direction . 


1 271 

l^(0;p,<p) = ^ yy' I.(0;-ia',cp)Sy(T;p,cp;-n'.cp')dp' dcp’ (282) 

0 0 


1 27T 

I^(T;-|i,cp) = yy* I.(0;-p',9)Ty (t; ~n,cp; -p' , 9 ')dp' dcp' (283) 
0 0 

1 2t: 

1^(0; p,cp) = ^ yy* I.(t;p, 9 )Tj^(t; p, 9 ;p',cp')dp' d 9 ' (284) 

0 0 


1 2 tz 

I^(0;-p,9) = yy* I.(t; p,9)Sj^(t; -p,9; n' ,9')dp' d9' (285) 

0 0 


These equations define the relationships of the bidirectional scattering func- 
tions S,, , T^t , St , and T,. for a layer i thick . The factor 1/ u is included so 
U U Li L 

that the scattering functions will be symmetric under interchange of directional 
angles. This operation could also be accomplished by multiplying the integrand 
by p' . Then, the equations would be identical to the customary defining equations 
for the bidirectional scattering (reflection or transmission) used for planar scatter- 
ers as found in the work of Love (ref, 4) and Nicodemus (refs. 5 and 6) . 

At first , the solution to the canopy problem seems to be found in solving for 
S^. This premise would be true if the canopy were a single layer; that is, if the 

canopy could be considered to be infinitely thick or if it were bounded on the lower 
side by a perfect absorber in the wavelength of interest. If radiation from any source 
is incident on the lower surface of the layer, then T,. must be known to find the con- 

tribution to upwelling intensity at t = 0 . If the upwelling intensity at the lower sur- 
face of the layer arises from reflection of the transmitted intensity of the layer, which 
is likely , then Tt-, and S,. must be known . Therefore , for multilayers , it is nec- 

essary to know the four scattering functions for each layer together with the rela- 
tionship that describes the attenuation of incident flux to give reduced incident flux. 
The set of equations that enable solving for these quantities follows . 
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Boundary Conditions 


Because the bidirectional scattering functions for the layer are properties 
of the layer itself, the functions are independent of the conditions of incident inten- 
sity and the bounding conditions outside the layer . Thus , the illumination and 
boundary conditions will be chosen in the most convenient way. 

Incident intensity is taken first as uniform, parallel, downward radiation in- 
cident on the upper surface . This incident radiation is represented by Dirac delta 
functions as 


1.(0; -p,<p) 


tiF 



(286) 


where is the magnitude of incident intensity . Incident intensity will then be 

taken as uniform, parallel radiation incident upward on the lower surface at t = t' . 
This relationship is also represented by Dirac delta functions as 


I. (t'; p,cp) 


tiF. 8 




(287) 


Using equations (286) and (287) in the defining equations for , T^, Sj^, and 
T gives 

Xj 


1(0; p,<p) 


F 

4il 



^.9; 


qi_,cp 


■) 


I(t';-p,cp) = ^ Ty ^T’;-p,tp;-p ,cp ^ 


and 


l(0;p,cp) = ^ 


F 

I(t’;-p,cp) = ^ Sj^^T’;-p,cp;p^,tp^^ 


(288) 


(289) 


(290) 


(291) 
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It is convenient to assume the layer is bounded by perfect absorbers or infi- 
nite space . Then , for either illumination 


l(0;-p,cp) = I(T';|j,,cp) = 0 


(292) 


where 0 ^ p < 1 . That is , there is no scattered , diffuse component incident on the 
layer either from above or from below , Note the existence of symmetry relationships 
between upper and lower surface transmission functions T^ and Tj^ . Thus , it 

is redundant to distinguish them; however, for simplicity within the derivation, 
the distinction will be maintained in the notation . 

Under the conditions of illumination specified, the radiation field at any level 
t within the layer can be separated into three parts: an upward scattered compo- 

nent I(t; p,cp) , a downward scattered component I(t; -p,cp ) , and an unscattered 
component composed of the reduced incident flux. The relationship describing 
the attenuation of or F_ is obtained from equations (276) and (278) with the 

omission of the term F. 

Integration is immediate to give exponential attenuation by a factor of e , 
where t is the absolute value of the optical thickness of material through which the 
incident flux passes . Thus , at any level , the reduced incident flux for F^ will be 
-(t’-t)/p_|^ -x/\i_ 

xF^e and the reduced incident flux for F_ will be xF e 

Parameterization Principles 

The basis for the technique employed to solve for the bidirectional scattering 
functions was developed by Chandrasekhar in solving for atmospheric transmission. 
This technique will be modified and expanded for use in canopy modeling. For 
downward radiation incident on the upper surface F_ , the following four principles 
can be formulated . 

1. For a layer t’ thick , the intensity I(x;|o.,cp) in the upward direction at 

-T/p_ 

any level results from the reflection of the reduced incident flux xF_e 

within the layer, and the reflection of the downward diffuse radiation I(t;-|j,, 9) 
(where 0 < p, < 1) incident on the surface t, by the layer of optical thickness 
t' - T, below T. The mathematical expression of this principle is 

. F- ~pT / \ 

I(t;p,9) = 4 jj e Sy/t’ - t; p,9; -p_ ,cp_J 

1 2x 

+ yy* ~ i:; p.tp; “P* -p' ,9')dp' dcp' (293) 

0 0 
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2. For a layer x' thick , the intensity I(t;-|o., 9 ) in the downward direction 
at any level t results from the transmission of the incident flux by the layer of opti- 
cal thickness t , above the surface at x , and the reflection at this same surface of 
the upward diffuse radiation l(t; p,cp) (where 0 < p < 1) incident on x from below. 
The mathematical expression of this principle is 


I(t; -p,9) 


4|I 


1 2tx 

+ 4 ^ J' f Sj^(t;-p,cp;p',cp')I(t;p',cp’)d^' ^9* (294) 

0 0 


3 . For a layer x' thick , the diffuse reflection of the downward incident radi- 
ation by the entire layer is equivalent to the reflection by that part of the layer of 
optical thickness t above the level t and the transmission of the diffuse intensity 
l(r, p,cp) incident on the surface x from below. The mathematical expression of 
this principle is 


F_ 


5u^T';p,9;-p_,9_^ 


F_ 

4p 




+ e I(t; p, 9 ) 


1 2tx 

^ 4i^ f f TL(t;p.9;i^'.9')I(x;p',9')d^' d9* 

0 0 


(295) 


4 . For a layer x' thick , the diffuse transmission of the incident intensity 
by the entire layer is equivalent to transmission of the reduced incident flux 
-x/p_ 

7xF_e and transmission of the diffuse radiation I(x;-p,9) (where 0 < p < 1) 

incident on surface x by the layer of optical thickness x' - x below x . The mathe- 
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matical expression of this principle is 


F_ / X F_ 

4fl 4jl ® - i;-p,cp;-p_,9_ 

t'-t 


+ e ^ I(t;-p,cp) 


+ 


1 

4n[i 


1 2n 


ffV'’ 

0 0 


t;-p,cp;-p' ,cp')I(t;-p' ,cp')dp' dcp' 


(296) 


The corollary set of four principles for radiation incident from below on 
surface t' can be formulated as follovvs. 

1. For a layer i' thick, the intensity I(T;-p.,cp) in the downward direction 
at any level t within the layer results from the reflection of the reduced incident 

flux TxF^e and the reflection of the upward diffuse intensity I(x;p, 9 ) 

(where 0 < p < 1) incident on the surface x, by the layer of optical thickness x, 
above x . The mathematical expression of this principle is 


I(x;-p, 9 ) = ^ e 


^x;-p,9;p^,9^^ 


1 2ti 

+ 4 ^ yy' Sj^(x;-p,9;p',9’)I(x;p',9')dp' d9* 
0 0 


(297) 


2. For a layer x' thick, the intensity I(x;p,9) in the upward direction at 
any level x results from the transmission of the incident flux by the layer of optical 
thickness x' - x below the surface at the level x , and the reflection at this same sur- 
face of the downward diffuse radiation I(x;-p, 9 ) (where 0 < p < 1) is incident 
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on t from above . The mathematical expression of this principle is 


I(T;^,cp) = Tj 


e 






1 2n 

+ y* y* Sy(T* - t; p,9; -p' ,cp')l(t; -p' ,9')dp’ d9’ (298) 

0 0 


3 . For a layer t' thick , the diffuse reflection of the upward incident radia- 
tion by the entire layer is equivalent to the reflection by that part of the layer of 
optical thickness t’ - t below the level t, and the transmission by this layer of 
the diffuse intensity I(t;-p,cp) (where 0 < p < 1) incident on the surface t from 
above . The mathematical expression of this principle is 


F F \ 


T'-T 

e I(t;-p, 9 ) 


1 2n 


^ yy'Tyft* - t;-p,9;-p',9’)l(T;-p’,9’)dp' d9* 

(299) 


0 0 


4. For a layer x' thick, the diffuse transmission of the upward incident in- 
tensity by the entire layer is equivalent to the transmission of the reduced incident 
-(t'-t)/p^ 

flux TiF^e and the transmission of the diffuse radiation I(x;p,9) (where 

0 < p < 1) incident on the surface t by the layer of optical thickness t above x. 
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The mathematical expression of this principle is 


x'-T 


4il ^ e ^x; li.cp; + e ^ I(x;p 


.<P) 


1 2tc 


4^ y * f d<p' 


0 0 


(300) 


System of Equations in Functions S and T 

The next step is the development of a system of equations for the bidirectional 
scattering functions . The technique employed was originated by Chandrasekhar 
and differs from section 51 of his "Radiative Transfer" (ref. 1) only in the addition 
of subscripts U and L and in the substitution of x' for . Equations deriving 

from principles expressed in equations (293) to (296) will be done simultaneously 
but separately from those derived from equations (297) to (300) . For convenience, 
set = -p_. Then, differentiate equations (293) to (299) with respect to x, apply 

the boundary conditions in equation (292), take the limit as x = 0 for equations (293), 

(296) , (298) , and (299) , and take the limit as x = x' for equations (294) , (295) , 

(297) , and (300) . For equations (293) to (296) , this procedure gives 


- ‘'-[l L , ) _ 

L dx 4p F_.<p_y 0^, 


1 2ti 


4 // d^. d<p' 

0 0 


(301) 



F- ^x’;-p,9;-p_,<p_J 

4p a? ^ 


+ 


1 

4icp 


1 2ii 

y * f Sj^(x'; -p,cp; p' ,cp') 
0 0 


dl(x;+p' ,cp') 
dx 


x=x' 


dp' dtp' 

(302) 
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and 


O-JiT 8f [ dr J,,,. 

1 271 p 

■"4^// dii' dcp' 

0 0 *- ^ - ^ 


0 = 7- 

4^ 


-jj-Tu^t';-n,<p;-^_,cp_^ 


aTy(T';-n,cp;-[i_,(p_] 


9t' 


-- r 

u dl(r,-^,y) 
^ ® dt 


1 27t 


T=0 


ik //Tu<t'-F.<p;-F'.<p') 


0 0 


dI(T; -|g' ,cp') 
dt 


t=0 


d|j,' dcp' 


Similarly, for equations (297) to (300), this procedure gives 




9SL(t';-ii,cp; iA+.fp+)' 


9t' 


1 27t p _ 

4 // 51M d,. 

0 0 L J ^ ^ 


dl (r, +p, 9 ) 

9Tfi 

dt 

t=0 ' ^ 


9t' 
1 2ti 


4^ If O'*' 


(303) 


(304) 


d9' 

(305) 


d9' 

(306) 
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and 


0 = 


^ L J t=0 

12 % p -| 

ik If "’ - d^ J <*'<’' 


(307) 


F 


+ e 




9TL(T';^,fp;^+>tp+) 
* Bd 


_■- 


dl(i; +|g,cp) 
dt 


1 2 % 

^ ^ 0 0 


X 


dl(t; +n' ,9') 
dt 


T=t' 


d|j.’ d9' 


(308) 


The derivatives in the preceding equations can be obtained from the previously 
derived equation of transfer and the expression for F in equations (276) and (277) . 
This procedure gives the set appropriate for equations (293) to (296) . 


dI(T; +[i,9) 
dr 


T=0 


= U 


^r';p,9;p_,9_^ 


4p Sur;ii>9;^_>9_) - F_(0 ;+p,9) 


(309) 


dl (r; -pi, 9) 
dt 


= F_(0;-^i,cp) 

t=0 ^ 


(310) 


dl(t; -^^,9) 
dt 


t=t' 


= F_(t’;+p,9) 


(311) 


dl (r; -p,9) 
dr 


T=r' 


4p 


?y^t';-p,9;p_,9^ - F_(r';-p,9)j (312) 
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where F_ is given by 


F_(0; |i,cp) 




,9; 





+ 


_L 

471 


1 2n 

J y' P(H,cp; 11 " ,cp")Sy /r'; ji” 

0 0 V 



and 


F_ (t'; ^i.cp) = j F 


T' 




,cp‘> 


-P_ 



1 2it 



The set appropriate for equations (297) to (300) is 


dl (t; -p,cp) 
dr 


T=T' 


i JL 

p 4p 




1^.9; 1^, 



F_|_ (t'; -p,9) 


dl ( t; -p,9) 
dt 


r=0 


+ ^ F^(0;-p,9) 


dl(t;+p,9) 

dt 


= ~n: 

t=t' ^ 


dl (t; +p,9) 
dt 


J t=0 




~F 

4jr Tl^^';P- 9;P+.9+^ - F+(0;p,9) 


(313) 


(314) 


(315) 


(316) 


(317) 


(318) 
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where is g^ven by 


F^(t';n,cp) = 


1 271 


^ PCp-j^p; l-i^ dcp’ 

H' 


(319) 


F,(0;^i,cp) = 


-I_ 
e p 


1 


(320) 


Combining equations (309) to (314) with equations (301) to (304) and com- 


U’ -"U’ 


/I + l\c /.I \ 

(5 5:jSu('‘^‘•9i-^»..cp^ + gp L 


= P 


|^,9;-|j,_,9_^ 

1 271 

4^ / ,9\^'dcp" 

0 0 \ / ^ 


1 271 


47t J'J' »‘P')p/~P' »«P'; ~^_.cp_') ^ d<p» 

0 0 V / ^ 

1 271 1 2ti: 


00 00 



as 


,cp j 


at' 


= p 


|jx,9;-ti_,cp ^exp 


T’ 1 2n 


^ ff V'.q>")Ty/i’;-p”,<p";-p_.q>_)^ dcp" 

n n ' / 


0 0 
t’ 1 2n 


^ JJ Tj^(T';p,cp;p’,cp’)p/p',cp';-p ,cp dcp' 
0 0 ' ' 


1 2n 1 2x 


;^//// Tj^(t'; |i,(p; p' ,cp’)P(p' ,cp';-p" ,cp") 

1 Ox i\ f\ n n 


0 0 0 0 


X Ty^'; p",cp"; p_>9_^^p~ dcp' cJcp" 


(322) 


i Ty^t';-p,cp;-p_,cp ^ 


^ dTufT';-p,cp;-p_,cp_^ 


at’ 


_ r; 

- P 


^-p,cp;-p ,9 ^ 

1 2x 

^ If P<-^^’‘P-’-J^^9'')Tu('t';-p''.9'';-p_,9_\^ d9" 

on V / ^ 


- 1 2x 


4x J' J' ' p,9» P* >9'> d9' 


1 2x 1 2x 


ibJJU Sj^(t'; -p,9; p' ,9')P(p' ,9'; -p" ,cp") 

1 Ox n f\ n rv 


0 0 0 0 


X Ty^t'; p",9"; p_,9^-i^ d9' d9” 


(323) 
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at' 


_ t; 

= e"^^ P/- 


^-^,cp;-|a_,9_^ 


-- 1 271 


ff P(-^,9;p".<P")Su/t';^",9";-^_,9_\^ 
n n V / 


0 0 
1 2tc 


y Jf Ty(t';-^,9;-p’ .9')P^-P' .9';-p_.9_'\^ 
0 0 ' 


1 27t 1 27t 


JJ JJ Ty(t'; -|X,9; -|i’ ,9')P(-n' ,9'; |i",( 

1 OTC n n n n 


0 0 0 0 


X Sy |t'; [i" , 9 "; ~n_ < 19 ' d 9 " 


3 S 


1 27T 


= + y yy' p(-p.9;-p'.9')SlA;-p’.9';p+.9+\^ 

^ '00 ' ' 


1 2t: 


+ ^ yy S^^(t';-^l,9;|a^9’)p/^l^9';^l^,9 d9* 

0 0 ^ ' 


1 271 


1 27T 


— ^ yy Sj^(t'; -^ 1 , 9 ; , 9 ’) yy P(n',9';-n",9") 

16 TI n n n n 


0 0 


0 0 


X s, 


^^,9^^^^ d9" ^ d9’ 


d9" 

d9' 

)") 

(324) 

• d9' 

(325) 
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0T 


T' 


= e 


^ P^^.9; ki+.9+^ 


- 1 271 


^ 4 ^ f f ,(p';[i^,(p\-^ d 9 ' 

on ^ / 


471 

J 

0 


1 

J_ 

47t 

/ 


0 ( 

1 



^ yy* Tj^(t';|i.cp;^',<p’)p/n',cp';n^,cp dtp' 

n n \ ' 


1 27T 1 27t 


(47t) 


b//// Tj^(t'; n.cp; n' , 9 ')P(n' , 9 ’; -^", 9 ") 


0 0 0 0 


X s 


L 

9T 


^t';-H”,9";|a^,9^^ d9* ^ d9" 


( 326 ) 


i ^.cp; 


= e 


JL 




1 27t 


4% 


J J Sy(T';^,9;-|i'.9')p/-n',9';n^,9 d9' 

n n V f 


0 0 
1 27t 


^ ff P(^A» 9 ;p^ 9 ’)TJ^^t';p^q>';^l+. 9 +^^ dtp' 
00 ' / 


1 271 1 27t 


(47t) 


bffff Sy (t’; n, 9 ; -ji’ , 9 ')P(-n’ , 9 '; ji" , 9 ") 


0 0 0 0 


X T, 


^t’;h",9";h^, 9^^^ d9" ^ d9* 


( 327 ) 
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= e 






1 271 


ff TL/t';n’,cp';^+.tp+V(-^.9;M'.<P’)^ ^9' 

0 0 ^ ' 


\i 1 2 ti 

• 1 - 


4tc 


J y* Ty(t';-n,<p;-ix',cp')P ^Ia'>9';^+.9+^^ ^9' 


0 0 
1 2tc 1 2ir 


( 471 ) 




0 0 0 0 


X Tj^^t';|a",cp";|^^,cp^^ dcp' dcp" 


(328) 


The derivatives can be eliminated by simple subtraction, which leaves four integral 
equations . The solution to these integral equations results in an expression for 
Sy . Using Sy , upwelling radiation for any specified condition of illumination can 

then be calculated by means of equation (282) . 
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TABLE I.- TOTAL REFLECTION AND TRANSMISSION FUNCTIONS 


p^(k,k')* 0 < ^1* < 1 0 < ji'* < 1 


Schematic 


Ty(k,k’)* -1 < ^1* < 0 0 < ^’* < 1 




p_ (k,k’)* -1 < < 0 -1 < ^’* < 0 

Li 


(k.k')* 0 < ^A* < 1 -1 < ^'* < 0 

Lt 




TABLE II." SLOPE AND AZIMUTH RANGES FOR PLANAR ELEMENTS 
CLASSIFIED BY INCIDENT OR VIEWING SURFACE 



Range of 
0. or 0 


0 < e. < I 


0 < < i 


Range of 0 


0 < 0 < ^ - 0. 
n 2 1 


? - ^ < 2 


Range of cp 


0 < o) < 2 tz 
^n 


(p. - 6. < cp < cp. + 6. 
1 ^n 1 


0 < e, < I 


i < % < t 


cp. + 6. < Cp < cp. + 27 t - 6. 
^1 1 ^n 1 


0 < ®s < 2 


0 < 9 < % - e 

n 2 s 


0 < cp <271 


0 < ®s < 2 


- e < e <7 

2 s n 2 


cp - 6 <q> <q) +6 

s ^n s 


0 < < 2 


^ - 9^ < 9^ < 
z s n 2 


cp + 6 <cp <cp +2 ti-6 
s s 


^ < 9. < 7C 


9. - ^ < 9 <5 

1 2 n 2 


cp. - Tt + 6. < cp < cp. + 7t - 6. 

1 ^n 1 


^ < 9. < X 


0 < 9„ < 9. - 7 
n 1 2 


0 < cp < 27t 


I < 9i < . 




Cp. + TT - 6. < Cp < Cp. + 71 + 6. 

1 ^n 1 




s 2 n 2 


<p. - Tt + < <p^ < CPg + ,1 - 63 




0 < 9„ < 9^ - 
n s 2 


0 < cp < 27t 




®s ■ I ^ ®«n I 


cp +71-6 <<P <<P +7C + 6 

s n s s 




I.- PHASE- FUNCTION INTEGRATION LIMITS 


’ll- 

“si 

'll- 

c 

(0 

to 

'i2- 

“si 

’i2- 

“s2 

’ll- 

“s 

'l2- 

I-s 


Bidirectional 

function 

involved 


0 < 0„ < 7 0 < 0^ < 


n ^ 


0 < cp < 2ic 


» < ®s < ? 5 - ®s < ®n < I ■ •Ps ■ ®s < ^ 'Ps * 


0 < e. < i i - 0 j < e_ < 


? I- e, r 


■ ®i < '^n < * ®i 


% < T - ®i «>s ^ »s < "Pn -Ps * - ®8 


5 I- ej ^ ^ 7 


< „ < hi * ®i 

■Ps * ®s " hs * ■ ®s. 


® < ®s < I I ' ®i < ®n < ^ ■ ®s >Pi * 6j < <p„ < <Pj * 2n - 6j 


“ ^ ^ 5 I ' [ej ^ ”n ^ 7 

0 < »s < I ? - Pj] < ®n < I 

I<e,<n e 3 - 5 <e„<|- 9 j 


<p + 6. ) Ftp. + 2n - 6."| 

hi * ®i < „ < hi ■ ®il 
K * ®s " Ts ■ ®sj 


" ®s ■ 1 < ®n < I ■ ®i 'f’s ■ ” * ®s < ‘I’n < fs ^ - ®s 


<®n<^ 


<0 <7l O<0^< 


•Pi ■ ®i 

<l>. - ” + 6 


0 < q> < 2 ti 


hi * ®i 

'Pn < [<Ps + n - 6, 


7 < 0s < " ®s - I < ®n < I ■ 0. <Ps * ” - ®s < “Pn < "Ps * " * ®s 


<e<s - ‘ <e„< 

e. - 5 " 


Pj_(k.k') 0<9j<^ 


O< 0^<7 O< 0 _< 


(p. - s. r<p. + 6 . 

■Ps * ” ■ ®s ^ ^ L<P3 - ” * 

<Pi * ®i < «, < hi * ■ 

<Ps ■ ” * ®s 0 hs * ’' ■ 
(p. + 5. < <p < cp. + 27t - 6. 


”n ^ “s ■ I *Pi ’ °i 'Pn ^ 'Pi 


<p. + 6j 

<Po + " ■ 


0 < cp < 2 tc 


[cp. + 2 ti - 
[tPg + ’I + 


Tj^Ck.k') J < 0 . < 7 t 

tj^Ck.K') |<0. <7t 


L"i ■ 

0 < < I I - < ®n < ®i - I ‘Ps - ^ < ^n < 

® < I ®i - I < ®n < ^ ■ ®s q>i ^ n - 6. < cp^ < cp. + 7. + 6. 


cp - 6 fcp + 6 

rs s < to < ® s 

|cpi + 7 t - 6 . '^n ^cp. + 7 t + 6 . 


P, (k.k') I ^ < 0. < 7t 0 < Og < I 7 - Qg < ®n < ®i " I + 6„ < cp^ < cp„ + 275 - 6^ 



TABLE III.- PHASE- FUNCTION INTEGRATION LIMITS - Concluded 
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0<u. <1 0<u<l lto 4 5, 6 9 7, 8 

0 < n- < 1 -1 < n < 0 10, 11 12 to 15 17, 18 16 

1 s 

-1 < n- < 0 -1 < n < 0 36 34, 35 28 to 31 32, 33 

1 s 

-1 < ji. < 0 0 < n < 1 25, 26 27 23, 24 19 to 22 

X s 

See table III, first column. 


TABLE V.- POLAR COORDINATES WITHIN THE 
INTEGRATION FRAME 


Region 

0 

9 

a 

0 < 0 < It 

7t < (p < 2it - 0 
^ n 

b 

0 < 0 < It 

% - Q < (p < Tt 

n ^ 

c 

0 < 0 < Tt 

2it - 0 < cp < 27t 

n ^ 

d 

0 < 0 < 7t 

0 < cp < Tt - 0 
^ n 
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